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LO SUMMARY 


A finite difference method for solving the unsteady transonic flow about harmonically oscillating wings is 
investigated. The procedure is based on separating the velocity potential into steady and unsteady parts and 
linearizing the resulting unsteady differential equation for small disturbances. The differential equation for the 
unsteady velocity potential is linear with spatially varying coefficients and with the time variable eliminated by 
assuming harmonic motion. 

The work of this report is a direct extension of earlier studies and includes correlation with experimental 
results for two rectangular wings and investigations of possible solution techniques for three-dimensional wings 
of more general planform. 

The main results of the study are as follows: 

1. An alternating direction implicit (ADI) procedure is investigated and a pilot program is developed for both 
two- and three-dimensional wings. This program provides a relatively efficient relaxation solution without 
previously encountered solution instability problems. 

2. Pressure distributions for two rectangular wings are calculated and the results correlated with experimental 
data. 

3. Conjugate gradient techniques are developed for the asymmetric, indefinite problem. The conjugate gradient 
procedure is evaluated for applications to the unsteady transonic problem. Several preconditioning methods 
are investigated. 

4. Difference equations for the alternating direction procedure are derived using a coordinate transformation 
for swept and tapered wing planforms. For a coordinate transformation which is continuous up to the second 
derivative, the ADI method converged for sweep angles up to 45 deg. The pressure distribution for the swept, 
untapered flat plate is correlated with the kernel function method. 

The results for the conjugate gradient method are preliminary and further testing of the method is contem- 
plated. Of the techniques studied, the most efficient procedure for the analysis of three-dimensional wings on 
computers of limited memory (i.e., on the order of one million words or less) appears to be the ADI method. 



2.0 INTRODUCTION 


The development of a capability for calculating unsteady transonic airforces for use in flutter analysis con- 
tinues to be of interest. The considerable effort on this subject falls generally into two categories; (1) finite dif- 
ference solutions of the transonic small perturbation equation for the velocity potential and (2) a modified strip 
theory procedure in which the subsonic coefficients are corrected empirically for transonic flow characteristics. 
In the first category, there are two main approaches; the first consists of a time integration of the nonlinear 
differential equations for the unsteady velocity potential, and the second involves separating the velocity poten- 
tial into steady and harmonically varying unsteady parts. The latter is the approach we have used, and its for- 
mulation results in the usual nonlinear differential equation for the steady velocity potential and a linear equa- 
tion for the unsteady velocity potential, with spatially varying coefficients which are functions of the steady 
velocity potential. Developments following this last approach have been documented in a series of NASA CRs 
(refs. 1 through 7). Results have shown that the two-dimensional, typical section problem can be handled by this 
harmonic approach relatively efficiently and accurately. However, the direct numerical solution technique 
which is successful for two-dimensions appears to be expensive when used on the full three-dimensional 
problem. 

The purpose of the work discussed in this report is to explore the feasibility of solving the three-dimen- 
sional problem. This was accomplished by investigating several different procedures. The first is a finite dif- 
ference formulation which results in a set of equations to be solved by ADI relaxation techniques. Then conju- 
gate gradient techniques were investigated. Finally, the out-of-core direct solution module was rewritten from 
two dimensions to three dimensions and transferred from the CYBER to the CRAY. 

The ADI method was derived by expressing the time-dependent differential equation for small perturba- 
tion transonic flow in difference equations using the alternating direction implicit (ADI) technique and then 
making the assumption of harmonic motion. The resulting ADI procedure does not have the frequency limita- 
tion on convergence, as do the standard block relaxation techniques. The ADI method was first tried on the two- 
dimensional problem. Correlations of examples for airfoils of vanishing thickness and for finite thickness airfoils 
are presented in section 5 with some discussion of convergence rates. The method was also extended to the rec- 
tangular wing for which the equations are derived in appendix A. Results for configurations of vanishing thick- 
ness are correlated with results from the integral equation method which is described in section 6 of reference 19. 
Also, results for rectangular wings with a 5%-thick circular arc airfoil and with a 12%-thick supercritical airfoil 
are correlated with experimental measurements. Problems in convergence were encountered when the method 
was applied to swept wings using a coordinate transformation aligning the streamwise variable with the lead- 
ing and trailing edges of the wing. This was alleviated by using a swept wing coordinate transformation having 
continuous second derivatives across boundaries of mapping regions. Solutions for the swept untapered wing 
were obtained for sweep angles up to 45 deg. and correlated with the integral equation method of references 8 
and 9. The solution for a 50 deg. swept untapered wing failed to converge. 

The classical conjugate gradient method is for solving systems of linear equations of the form Ax = b with 
coefficient matrices which are symmetric and positive definite. This procedure proved to be efficient for large, 
sparse, well conditioned systems. Work for this report has resulted in an algorithm for system matrices which 
need only to be nonsingular. Before applying the conjugate gradient algorithm, the coefficient matrix for the 
transonic problem must be preconditioned in order to speed convergence. Three methods of preconditioning were 
tried and are discussed in section 7. The most effective procedure is derived from the partial LU decomposition of 
the coefficient matrix. The derivation of the complete algorithm is presented in appendix E. 
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The pilot program for the out-of-core direct solution of two-dimensional airfoils was extended to rec- 
tangular three-dimensional wings. This pilot program was also converted from the CYBER to the CRAY com- 
puter to take advantage of the significantly larger core storage. The resulting test runs were expensive, but the 
cost was dominated by charges for input/output operations. The availability of machines with even larger core 
storage capability, and the characteristic that a large number of mode shapes may be handled for a minimal 
additional cost over that for one mode, warrant further work on the direct solution. 
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3.0 ABBREVIATIONS AND SYMBOLS 

a Amplitude of wing oscillations 

b Root semichord 

Cp Pressure coefficient, (p-po)/(^/^ PoUq^) where p is the local pressure, Po the freestream static pressure, and 

Po the freestream air density 

A Matrix of coefficients, also aspect ratio (Section 6) 

f Frequency in Hertz 

fo Undisturbed wing or airfoil shape 

fi Unsteady contribution to wing or airfoil shape 

Y x,y subscripts and indices for points in the mesh 

i VT ^ 

I Identity matrix 

k Reduced frequency based on semichord, 2n fb/U; same as co. 

K Transonic parameter, (l-M^VCM^e) 

le Leading edge 

M Freestream Mach number 

n (n^, Uy, r^.), unit normal vector to shock 
q co^/e - ico(y - l)(Po 

XX 

t Time in units of b/U 

te Trailing edge 

u K-(y + 1)(Po 

X 

U.Uq.V Freestream velocity 

X Freestream coordinate. Vector of unknowns in matrix equations. 

Xo,yo Physical coordinates, made dimensionless with the root semichord 
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x,y Scaled coordinates (xo,HYo) for the two-dimensional problem 
Xq Steady chordwise shock location 

Complex amplitude of shock oscillation 
Magnitude of X^ 
a Angle of attack 

P 2 = <u^> 

Y Ratio of specific heats for air 

ACp Jump in pressure coefficient across airfoil or wake 

At Time step in ADI procedure 

A<pj Jump in (p^ at plane of wing or vortex wake 

A(Pi Jump in (p^ at wing trailing edge 

te 

8 Thickness ratio; also finite difference operator 

e (8/M)2^3 

toM/(l-M2) 

p Scale factor of Yq, p = 8^3 ^2/3 

rf Fraction of semichord 

Swept wing coordinates 

(}) Unsteady time dependent perturbation potential 

(p , (p Complete, scaled perturbation velocity potential; also used for the unsteady potential in finite difference 

equations with subscripts 

(pQ, ifQ Steady scaled perturbation velocity potential 
(Pi , <pi Unsteady scaled perturbation velocity potential 
0) Reduced frequency in radians; same as k 
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Notation 


[ ] Denotesjump in quantity across shock 
< > Denotes mean value of quantity at shock 
A Denotesjump in quantity across airfoil except for At 
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4.0 FORMULATION AND SOLUTION 


Since the mathematical derivation of the method for the solution of the unsteady velocity potential for the 
flow about a harmonically oscillating wing is presented in reference 1, the discussion here will be limited to a 
brief outline of the procedure for two dimensions. The complete nonlinear differential equation was simplified by 
assuming the flow to be a small perturbation from a uniform stream near the speed of sound. The resulting 
equation for unsteady flow is 

[K - (7 - 1) - (7 + 1) v^x] V’xx ^ ‘^yy - (2<^xt ^t) ° (4-D 

where K = / QVI^e), M is the freestream Mach number of velocity in the x-direction, x and y are made 

dimensionless to the semichord b of the airfoil and the time t to the ratio b/U^. With the airfoil shape as a function 
of time defined by the relation 

yg = 5f(x,t) 

the linearized boundary condition becomes 


iPy = fx (x,t) + (x,t) (4.2) 

The quantity 5 is associated with properties of the airfoil (such as maximum thickness ratio, camber, or 
maximum angle of attack) and is assumed to be small. The coordinate y is scaled to the dimensionless physical 
coordinate y^ according to 

y = M^^^yg 

and 8 is given in terms of 8 by 

e = (6/M)2/3 

The pressure coefficient is found from the relation 

^p = ■ 2« (‘^’x + n) 

The preceding differential equation is simplified by assuming harmonic motion and by assuming the 
velocity potential to be separable into a steady-state potential and a potential representing the unsteady effects. 
We write for a perturbation velocity potential 

= <^0 (x.y) + v’l (x,y) 6^“^* (4.3) 


and for the body shape 

yO = 6f(x,t) = 5[fo(x) + fi(x)ei‘"‘] 

Since the steady-state terms must satisfy the boundary conditions and the differential equation in the 
absence of oscillations, we obtain 


[K - (7 + 1) = 0 (4.4) 

with 

= y = -1<X<1 (4.5) 
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On the assumption that the oscillations are small and products of cpi may be neglected, equations (41) and 
(4.2) with the aid of equations (4.4) and (4.5) yield 

{[K - (7 + 1) v’Ox] ^ ‘^lyy “ (2iw/e) + qv^i = 0 (4.6) 

where 

q = coV € - iw (7 - 1) v^Oxx 


subject to the wing boundary conditions 

" ^Iv -l^x<l (4.7) 

y ^ 

The boundary condition that the pressure be continuous across the wake from the trailing edge was found 
in terms of the jump in potential to be 


A<Pl = (^"^te) (4.8) 

where is the jump in the potential at x = just downstream of the trailing edge and is determined to 
satisfy the Kutta condition that the jump in pressure vanish at the trailing edge. The quantity Acp^ is also used in 
the difference formulation for the derivative cply^ to satisfy continuity of normal flow across the trailing edge 
wake. 

For the set of difference equations to be determinate, the boundary conditions on the outer edges of the 
mesh must be specified. In the original unsteady formulation, these boundary conditions were derived from 
asymptotic integral relations in a manner parallel to that used by Klunker (ref. 10) for steady flow. A later for- 
mulation (ref 3) applies an outgoing plane wave boundary condition to the outer edges of the mesh. This bound- 
ary condition is numerically simpler to apply and is equivalent to the first order nonreflecting boundary condi- 
tions derived by Engquist and Majda (ref. 11). 

A computer program for solving the steady state transonic flow about lifting airfoils based on equations 
(4.4) and (4.6) was developed by Cole, Murman, and Krupp (refs. 12 and 13). Steady-state solutions required for the 
coefficients of equation (4.6) were obtained by using the latest version of this program, TSFOIL (reference 14), 
and also by using the steady-state solutions from the time dependent method of Rizzetta and Chin (ref 15) called 
EXTRAN2. An additional improvement on locating the shock is also available both in TSFOIL and EXTRAN2. 
For airfoils at high Mach numbers and angles of attack, TSFOIL and EXTRAN2 predict shock positions consid- 
erably aft of experimental measurements. To overcome this difficulty, Jou and Murman (ref 16) developed a 
phenomenological model for the displacement thickness effects of shock wave boundary layer interactions. A 
wedge (or ramp) is introduced behind the shock to simulate the thickening of the boundary layer. This procedure 
was not needed for the results reported here, however. For three dimensions, we used XTRAN3S, the time 
dependent method of Borland, Rizzetta, and Yoshihara in reference 17, to compute the required steady state 
potential. 

The similarity of the unsteady differential equation to the steady state equation suggests that the method 
of Cole, Murman, and Krupp should be an effective way to solve equation (4.6) for the unsteady potential (Pj. Note 
that equation (4.6) is of mixed type, being elliptic or hyperbolic whenever equation (4.4) is elliptic or hyperbolic. 
Central differencing was used at all points for the y derivative and all subsonic or elliptic points for the x deriva- 
tives. Backward (or upstream) differences were used for the x derivatives at all hyperbolic points. The preferred 
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numerical approach to solving the resulting large-order set of difference equations is a relaxation procedure, 
which permits the calculation to be made as a sequence of relatively small problems. However, as discussed in 
preceding NASA reports by the authors (refs. 3 and 4), a significant problem of solution convergence with the 
relaxation procedure was encountered that severely limits the range of Mach number and reduced frequency for 
which solutions may be obtained. Accordingly, an out-of-core solver (ref. 18) was developed to solve the complete 
set of difference equations simultaneously, which for two-dimensional flow is relatively efficient. 

The size of practical three-dimensional problems is such that the out-of-core direct solver would appear to 
be very expensive. This and other solution procedures are under investigation. First, conjugate gradient tech- 
niques have been examined as a more efficient means for obtaining a direct solution. Since the matrix of coeffi- 
cients is neither symmetric, nor always positive definite, nor well conditioned, special procedures must be used 
in applying the conjugate gradient technique. Discussion of the algorithms tried during the current work is 
presented in section 7. Second, a relaxation-type method based on the time-dependent ADI method to obtain 
frequency domain solutions was derived which does not have the frequency limitations of classical block relaxa- 
tion. If, in place of substituting equation (4.3) into equation (4.1), we use instead 

v>(x,y,t) = <^o (x,y) + <Pi (x,y,t) 
we obtain the following linear differential equation for (p^: 

([K - (7 + 1) V>0x] ‘^Ox)x *^lyy - (2<^lxt + ‘^Itt) 

The ADI method of Borland, Rizzetta, and Yoshihara (ref 17) was applied to the equation, and then the harmonic 
assumption for the nth approximation 

j (x,y,t) = </> j (x,y) 

was substituted into the difference equations. Here At is the time step between successive approximations, and 
has no physical significance to the problem. This yields a set of equations which has the same form as alternating 
row and column relaxation but for which the solution is convergent for all frequencies. The finite difference 
procedure was derived in reference 7 for two-dimensional flow. The equations for three-dimensional flow are also 
derived in Appendix A for the rectangular wing and in Appendix B for the swept and/or tapered wings. 
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5.0 ADI RELAXATION SOLUTIONS FOR AIRFOILS 


In NASA CR 3537 (ref. 7), we showed that the direct solution of the difference equations for the harmonic 
motion of an airfoil in transonic small perturbation flow yields aerodynamic forces that take into account the 
effect of moving shock waves, at least in a first order sense. A program for two-dimensional flow based on the 
direct solution is reliable, efficient, and yields results which are in good agreement with experimental measure- 
ments (e.g., see ref 7). Also, typical section flutter boundaries calculated with airforces from this program exhibit 
the characteristic 'transonic bucket.’ The direct solution is also applicable to three-dimensional flow, but here 
the computer requirements are large and costly. Accordingly, a method of relaxation was developed in which 
simple harmonic motion was assumed after the finite difference formulation rather than before (see section 4.0). 
The resulting set of difference equations is solved by implicit alternating direction techniques, and hence the 
procedure has been entitled the “ADI” method. This ADI method does not have the convergence problem of the 
conventional block relaxation procedures. The differential equations were derived in NASA CR 3537 and pre- 
liminary two-dimensional examples were included. The derivations for three-dimensional rectangular wings 
and for swept, tapered wings are given in appendices A and B of this report. 

The equations of the ADI procedure and direct solution differ in several respects. First, there is a time step, 
At, which appears in the ADI equations but not in the direct solution. This time step, which has no physical 
significance, does lead to truncation errors which are discussed in the next section. Second, the shock point oper- 
ator which proved to be satisfactory for the direct solution had to be modified for use in the ADI difference equa- 
tions. This is discussed in section 5.2. 

To test the practicality of the method for eventual application to three dimensions, the two-dimensional 
ADI method was applied to a number of problems which include airfoil configurations both with and without 
thickness. Correlations are made with solutions from a program based on a compressible kernel function (see ref 
19), and from the direct solution program evaluated in reference 7. 

Because of the low cost of each two-dimensional iteration, the solutions were generally run until the max- 
imum of the difference between time steps was less than 10'^, a criterion found acceptable for the earlier relaxa- 
tion techniques. The number of iterations depends upon the time step chosen, smaller time steps requiring a 
proportional increase in the number of iterations. There also is a considerable Mach number influence on the 
rate of convergence. For a frequency of about 0.8 and choice of time step of 0.4, it was found that for M = 0.9 
convergence required about 600 iterations while for M = 0.7 less than 300 were needed, with the variation with 
Mach number being roughly linear. 

The following table gives data on the number of iterations required for convergence for several time steps 
for the severe case of M = 0.9 and a reduced frequency of k = 0.9 for a symmetric flow using a 72 by 60 mesh (2160 
mesh points). 


Time step 

Iterations 

CP seconds (CRAY) 

0.1 

1000 

55 

0.47 

26 

40 

0.6 

514 

28 

0.8 

406 

22 

1.0 

372 

20 
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An indication of the superiority of the direct solutions for two dimensions is seen from the time required for the 
full (72 by 60) and half (72 by 30) matrix solutions; that is, 17 cps for the full matrix and 4.5 cps for the half matrix. 

5.1 CORRELATION OF EXAMPLES FOR AIRFOILS OF VANISHING THICKNESS 

The ADI method was first applied to calculating the pressure distribution on a flat plate pitching harmon- 
ically about the leading edge. Agreement between solutions from the ADI method and the program of reference 

5.1 using an integral equation method were good for a Mach number of 0.9 and reduced frequencies up to 0.3 as 
seen in figure 1. However, agreement for a reduced frequency of 0.45 is not as good (fig. 2), and for a reduced 
frequency of 0.6 the agreement was found to be poor (see fig. 3), Reducing the time step improved the solution 
somewhat. 

Basically, the time step is irrelevant to the problem but it does lead to truncation errors. By adding equa- 
tions (A-5) and (A-6) to (A-7) and equating the three successive values of (p, that is, setting (p^ = (p^ = <p*^ + we 
obtain a difference equation which we can recognize as resulting in the following differential equation (see equa- 
tions (A-62) and (A-64)): 



where Pj = e * When the time step At goes to zero with u = K, this equation reduces in the limit to the 
classical linearized equation for harmonic unsteady flow. To check whether the program contained basic errors, 
this differential equation (without the (p^^ term) was solved by the direct method used for the classical equation. 
The resulting solution agrees very closely with the ADI solution as seen in figure 4. This verifies that the 
algorithm was properly programmed. A second order (p^^ difference which required saving an additional time 
step was also tried but introduced no significant improvement. 

Since the direct solution of the regular equations used a second order central difference for the (p^ term, a 
second order backward difference was derived and introduced into the ADI method. The derivation of this (p^ 
operator as well as the cp^^ operator enabling a change in the time step during a run are presented in Appendix D. 
Considerable improvement resulted, and now the ADI is in good agreement with the kernel function results. 
This is shown in figure 5 for a reduced frequency of 0.6. Therefore, to obtain sufficient accuracy at the higher 
frequencies we must use a second order difference for the first derivative in x. 

The effect of varying the size of the time step, At, for M = 0.9 and k = 0.6 is shown in figures 6 and 7, with the 
results plotted with an expanded ordinate scale. Figure 6 presents solutions from the program of reference 19 
and ADI results for the smallest and largest values of At tried. For the range studied. At = 0.1 to 1.5, the best 
correlation is obtained with the smallest value of At. Generally, the ADI calculations move towards the integral 
equation results in monotonic fashion as At is decreased. A reasonable compromise with accuracy for the sake of 
economy of calculations is found by using At = 0.6. 

5.2 CORRELATION OF EXAMPLES FOR AIRFOILS WITH THICKNESS 

A second set of examples was prepared for airfoils of finite thickness. The pressure coefficient distributions 
for the NAC A 64 AOlO airfoil, oscillating in pitch about the quarter chord, tested by NASA Ames, were calculated 
with both the ADI and direct solution procedures. Calculations were made for a Mach number of M = 0.85 and 
reduced frequency of k = 0.15 and for M = 0.86 and k = 0.4. The resulting pressure distributions are presented in 
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figures 8 and 9. Agreement is quite good in both cases except for the real parts in the vicinity of the shock. In 
figure 8, the real part of the pressure aft of the shock from the ADI program is significantly larger than that from 
the direct solution. In figure 9, the ADI results underpredict the pulse while overpredicting the pressures aft of 
the shock. The difference in results between the ADI and direct solutions may well be a difference in shock point 
operators. Since use of the shock point operator which is the exact equivalent to that in the direct solution causes 
divergence in the ADI solution, the shock point operator was applied only to the second derivative with respect to 
X in the ADI program and not to both the first and second derivatives as in the direct solution. 

Pressure distributions for a Mach number of 0.85 and for reduced frequencies of 0.25 and 0.4 are also pre- 
sented in figures 10 and 11. The distributions are presented for two ADI solutions and the direct solution. The ADI 
solutions differ in the representation of the first derivative with respect to x, one using a first order finite dif- 
ference and the other a second order. It is seen that there is relatively little difference between the two. However, 
in this particular case, the shock pulse of the ADI solutions differs considerably from that obtained by the direct 
solution. 

Since a singularity occurs across the shock in the pressure distribution, a more accurate method of deter- 
mining the influence of the shock point operator is to plot the difference in velocity potential across the airfoil. 
This jump is plotted in figure 12 for k = 0.25. Two ADI solutions are compared with the direct solution and are 
indicated by the legend on the figure. Again, one solution uses a first order difference for the first derivative with 
respect to x and the other uses a second order difference. For a shock of this strength, a distinct jump in the 
unsteady potential difference is observed, which can be related to the amplitude of shock motion by the shock 
relations derived in appendix A of reference 7. We note that there is very little jump in the real part of the 
potential difference at the shock for either ADI method compared with the direct solution. 

To obtain a better shock representation, shock boundary conditions were derived for application to the ADI 
method of relaxation, and the derivation is presented in appendix D. The expressions are similar to those 
obtained for the direct solution and described in NASA CR 3537. An example of the captured shock on the NACA 
64A010 at the Mach number of 0.85 is shown in figure 13. When the shock strength is defined as 

H=«i+lj - uy 

some improvement in the jump in the real part of the unsteady potential for the airfoil oscillating in pitch about 
the leading edge is obtained, but results for the imaginary part are worse than those from the shock point oper- 
ator (fig. 14). When the shock strength is defined by the more realistic value of 

the jump in the real part of the unsteady potential difference is in close agreement with the direct solution, and 
the jump in the imaginary part is somewhat less than that for the direct solution. The corresponding pressure 
distributions are shown in figure 15. In the direct solution, shock boundary conditions were found to lead to 
somewhat different results for the jump in the unsteady potential from those obtained by the shock point oper- 
ator. In reference 7, it was shown that the shock point operator for the direct solution lacked one term in the shock 
jump conditions. The application of the shock boundary conditions appears to be an acceptable method of repre- 
senting the effect of the shock in unsteady flow in the ADI relaxation method. 
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To study the effect of first- and second-order difference representations of the (p^ term on airfoils with thick- 
ness, the preceding configuration was run at M = 0.85 and reduced frequencies of 0.25 and 0.40. The results of 
the ADI method for both representations of (p^ are compared with the direct solution in figures 10 and 11. The 
greatest deviation occurs in the vicinity of the shock where the pressvure pulse is significantly underestimated by 
the ADI calculations. Inclusion of the second-order difference does not appear to be as important for the airfoil as 
for the configurations of vanishing thickness. 


13 



6.0 ADI RELAXATION SOLUTIONS FOR 
THREE-DIMENSIONAL WINGS 


6.1 CORRELATION OF EXAMPLES FOR 
RECTANGULAR WINGS OF VANISHING THICKNESS 

Pilot programs for three-dimensional unsteady transonic flow using the direct solution procedure 
(OPTRAN3) and the ADI relaxation method (OPTRAD3) have been applied to a zero thickness rectangular wing 
of aspect ratio 3 which corresponds to the planform geometry of the wing model of NASA TN D-344. Typical 
results from these two programs are compared with corresponding calculations from the RHOIV program of 
reference 8 in figures 16 through 18. These results are for a Mach number of 0.9 and a reduced frequency, based 
on the semichord, of 0.10, with the wing oscillating in pitch about the leading edge. The correlation of the results 
from the two finite difference procedures with the results from RHOIV is considered very good, even close to the 
tip. 


For the real part, results from the three theories match nearly exactly for the inboard two chords, and 
RHOIV lies slightly below OPTRAD3 and OPTRAN3 for the tip chord. For the imaginary part, OPTRAD3 and 
RHOIV match closely over the inboard two chords with the OPTRAN3 curve lying slightly below. For the tip 
chord, OPTRAD3 and RHOIV match closely over the front part of the smface (OPTRAN3 again lies slightly 
below), while OPTRAN3 and RHOIV match closely aft of the midchord (OPTRAD3 lies slightly above). The 
better agreement of OPTRAD3 and RHOIV is surprising since truncation errors due to the time steps At were 
observed in the two-dimensional solutions. 

The solutions are for a 50 by 20 by 40 xyz mesh (symmetry in z assumed so the order of the coefficient 
matrix is 20,000). As expected, the direct solution proved to be very expensive to run for this three-dimensional 
problem, and the ADI solution requires about a quarter of the computing resources of the direct solution. 

6.2 CORRELATION OF EXAMPLES FOR RECTANGULAR WINGS OF FINITE THICKNESS 

Steady-state pressure distributions for the wing of NASA TN D-344 with a 5%-thick circular arc airfoil are 
presented in figure 19. Included in the figures are pressure distributions for the upper and lower surfaces (as 
digitized from the graphs of TN D-344, ref. 20) and analytical results from XTRAN3S. It is noted that, despite a 
symmetric airfoil shape and zero angle of attack, the pressures for the upper and lower surfaces do not coincide. 
The report mentions a cross tunnel variation in stream angle which may account for the differences. The pres- 
sure distributions are presented for four spanwise stations. The analytical stations are closely matched to the 
experimental stations in all cases, so that a good estimate of the correlation between theory and experiment may 
be obtained by reviewing the figures. At the root (fig. 19), the analytical results extend to a larger negative pres- 
sure coefficient than the experimental results. Also, the analytical distribution shows a relatively sharp shock 
while the experimental results show some supersonic flow but almost no shock. Since there is a shock at mid 
semispan, a shock at the root would be expected too. The lack of a shock at the root may be due to the cross tunnel 
variation in stream angle and may also be due to boundary layer effects. The pressure orifices at the root are on 
the tunnel wall, rather than on the wing, and there is no mention of a splitter plate to remove the boundary layer 
on the tunnel wall. A thick boundary layer over the orifices might mask shock effects in the pressure 
measurements. 

At the mid semispan section, y/Ab = 0.5 (fig. 19), the correlation between theory and experiment is quite 
good. The analytical result matches that for the upper surface well, particularly with respect to the shock 
strength. The maximum pressure is slightly less for the analytical result than for the experimental result, other- 
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wise the analytical result seems to lie between the measured distribution for the upper and lower surfaces. For 
y/Ab = 0.70 in figure 19, the analytical results fall slightly below the measured distributions, however, the shock 
strengths appear about the same. For y/Ab = 0.9 in figure 19, analytical results fall slightly farther below the 
measured distributions. Again, the shock strength looks about the same for analytical as for the measured 
results. 

Using the potential distribution from the steady state solutions of XTRAN3S just discussed, we analyzed 
the wing of NASA TN D-344 (ref 20), with its circular arc airfoil section, under the bending mode presented in 
figure 20. The results are for M = 0.9 and reduced frequency of 0.13. Three distributions are presented in each 
figure: (1) the ADI relaxation calculations for the wing with the circular arc airfoil, (2) the ADI calculations for a 
flat plate and (3) the experimental measurements as obtained from figure 9c of TN D-344 using an automated 
digitizing process on the unflagged data. The amplitude and phase angle of the pressure difference across the 
wing are presented for three spanwise stations: y/Ab = 0.5, 0.7, and 0.9. 

The results from the ADI method in figures 21 and 22 for the thickness case reflect a much sharper and 
stronger pressure pulse than do the experimental measurements. Comparison of the amplitudes at the three 
spanwise stations is presented in figure 21. At y/Ab = 0.5 and 0.7, the experimental and anal}d:ical shock loca- 
tions agree relatively well. At y/Ab = 0.9, ADI calculations show a very small shock pulse, while the presence of a 
shock in the experimental data is unclear due to variations between the flagged and unflagged data in figure 9c 
of TN D-344. Ahead of the shock, the amplitude of measured pressure is greater than the calculated pressure for 
all three spanwise stations. 

A comparison of the phase angle distributions is given in figure 22. Overall, the correlation appears about 
like that for the pressure amplitudes. The analytical results show a sharp rise in phase angle across the shock at 
all three spanwise stations, with the phase angle decreasing aft of the shock for the two inboard stations. 
Although there are only two data points aft of the shock, the experimental data does generally resemble the 
analytic data. At y/Ab = 0.5 and 0.7 (fig. 22), the experimental data show a much larger jump in phase angle 
across the shock than the anal 3 d:ic data. Behavior of the experimental data aft of the shock, although described 
by just two data points, generally matches the analytic behavior at y/Ab = 0.9 in figure 22. The experimental 
data does show a large spike in phase angle (presumably due to a shock), while the analytic results show a jump 
across the shock. 

The results of analyzing the wing of NASA TN D-344 presented in figures 21 and 22 are replotted in figures 
23 and 24 to include the second set of experimental data (the flagged data) from the report. The flagged data were 
included in the report to show the repeatability of the measurements. Measurements in the vicinity of the shock 
show large variations between the two sets of data. 

Three-dimensional plots of the real and imaginary parts of the calculated pressure distribution for this 
problem are shown in figures 25 and 26. Figure 25 is for the vanishing thickness configuration, while figure 26 is 
for the wing with a circular arc airfoil section. The pressure pulse due to the presence of the shock shows up 
dramatically for the circular arc configuration. There is a noticeable blip in the pressure distributions at the 
trailing edge of the tip chord for both the zero thickness and airfoil distributions. Such a blip has been encoun- 
tered in other time-integration finite difference solutions. It seems to be associated with the solution mesh and is 
considered to make an insignificant contribution to the overall flow solution. 

The overall correlation between theory and experiment is somewhat less than satisfying, although the 
basic properties of the pressure distribution in the experiments are reflected in the theoretical calculations. The 
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theoretical solution has been checked for convergence by running an additional 100 iterations with no apparent 
change in pressure distribution. 

More recent experiments reported by Rickets, Sanford, Seidel, and Watson of the NASA Langley Research 
Center (ref. 21) correlate better with calculations by the frequency domain method. A rectangular wing of aspect 
ratio 4.0 with a supercritical 12%-thick airfoil is oscillated in pitch about the leading edge. Freon gas was used in 
the wind tunnel for which the ratio of specific heats is y = 1.131. 

The steady state solution for a Mach number of 0.7 and angle of attack of 2 deg. was computed by 
XTRAN3S. The resulting chordwise pressure coefficient distributions are shown in figures 27 through 30 at 
fractions of the semispan of q = 0.31, 0.59, 0.81, and 0.95, respectively. The agreement of the experimental meas- 
urements with theory is seen to be quite good. 

Using the steady-state potential just described for calculating the coefficients of the three-dimensional ADI 
method, we ran the program for the rectangular wing oscillating in pitch about the leading edge at reduced 
frequencies of 0.178 and 0.356. Figure 31 presents the distribution of amplitude of the jump in pressure across the 
wing for k = 0.178 at spanwise stations of Tf = 0.31, 0.59, 0.81, and 0.95, while figure 32 shows the corresponding 
phase angle distributions. The pressure amplitude over the aft 3/4 chord of the wing is in good agreement with 
experimental values. The OPTRAD3 results show the leading edge singularity at all four spanwise stations, 
while this singularity appears only for the outer two chords in the experimental data. A shock pulse appears just 
aft of the leading edge for the two inboard chords in both the OPTRAD3 and the experimental results. The shock 
pulse is sharper in the theoretical calculations, and the theory shows a singularity at the leading edge which has 
the same sign at all cross sections. This is not seen in the experimental values. The agreement of the theory with 
the phase angle measurements is not as good as for the amplitude. The experimental results indicate a rise in 
phase angle to about 3/4 chord and then a fall to zero, while the theory continues upward. Similarly, comparison 
of the ADI results with the first harmonic from XTRAN3S, as digitized from the graphs of reference 21, is shown 
in figures 33 and 34. The pulse near the leading edge from XTRAN3S is stronger than the harmonic solution for 
rf = 0.31. For the more outboard spanwise locations the agreement is better. The phase angle from XTRAN3S 
shows a rapid rise at the trailing edge not observed either in the experiments or the ADI solution. For k = 0.356, 
the correlation of the ADI solution with experimental results is similar to that for k = 0.178, as seen in figures 35 
and 36. The correlation with the first harmonic of XTRAN3S is shown in figures 37 and 38, and is similar to the 
results found for the correlation of two-dimensional results with XTRAN2. 

6.3 SOLUTION CONVERGENCE PROBLEMS FOR SWEPT WINGS 

The ADI program was modified to include the swept wing coordinate transformation. The problem of con- 
struction of a mesh system for a swept wing is considerably simplified by using a coordinate transformation 
which defines the leading and trailing edges as single values of the streamwise coordinate system. This makes it 
easier to refine the grid in the neighborhood of the leading and trailing edges. 

Since XTRAN3S was selected to compute the steady-state potential distribution required in the coeffi- 
cients of the difference equations, the mapping used in XTRAN3S was adapted to our program. The coordinate 
transformation is given by 


V = y 

^ = z (6.31) 
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where S(y) is the semichord at the span location y. This differs from the transformation in reference 17 in order to 
make the leading and trailing edges be given by ^ = -1 and 1, respectively. The differential equation resulting 
from applying this transformation and the resulting difference equations are derived in appendix B. 

The three-dimensional ADI relaxation program using the swept wing coordinate transformation has been 
applied to the clipped delta planform (as a zero thickness configuration) of Hess, Wynne, and Cazier (ref 22), 
which has a swept leading edge of 50.45 deg. an unswept trailing edge, and a semispan of 0.708 chords. The time 
step had to be reduced to 0.002 before convergence occurred. However, after 570 iterations, the solution again 
started to diverge. A plot of the pressure distribution at the time when convergence stopped was in very poor 
agreement with the RHOIV solution. A modified clipped delta wing with a leading edge sweep angle of 20 deg. 
was also run to determine the effect of the leading edge sweep. The convergence rate was somewhat improved. 

The grid region for the 50.45 deg. clipped delta wing is shown in figure 39. Note that the outboard section 
covers only a small region of the flow and may account in part, at least, for the lack of agreement with the RHOIV 
solution as well as the tendency toward solution divergence. 

To find the effect of sweep alone, a simple swept untapered wing was analyzed. Convergence was obtained 
for a sweep angle of 20 deg. but failed at 30 deg. When the mapping was changed to make the leading and 
trailing edges turn normal to the line of symmetry as in figure 40, the convergence was extended to 30 deg. but a 
40 deg. sweep failed to converge with the maximum growth occurring at the point where the edges curved 
sharply. At this point the mapping has discontinuous second derivatives. This suggests that a mapping with a 
continuous second derivative may improve the convergence. To test this, a swept untapered wing was analyzed 
using a cubic to bend the leading and trailing edges normal to the plane of symmetry. A finer grid spacing in the 
streamwise direction was also applied on the downstream boundary. It was found that, for the new transforma- 
tion, convergence was obtained for sweep angles up to 45 degrees. 

Additional studies concerning solution divergence of swept configurations using coordinate transforma- 
tions may be found in references 23 and 24. In the former, care is taken to assure that the derivative which 
appears as a coefficient in the transformed differential equation, is well behaved in the solution region. In refer- 
ence 24, linear stretching is applied to the transformation in the regions both upstream and downstream of the 
wing planform to assure that the region of perturbed flow is included in the flow solution region. Both concepts 
appear to significantly improve convergence. We assume that our new transformation is close to that used in 
reference 23. However, we have also made a point of removing the singularities resulting from the sharp plan- 
form apex at the root by artificially rounding off the leading and trailing edges at the root. 

To test convergence at the higher frequencies, the 45 deg. sweep wing was calculated using k = 0.5. The 
solution converged as easily as for the lower frequency calculations. As seen from figures 50 through 52, the 
agreement with RHOIV is, however, not as good. As in the two dimensional solutions, the second order difference 
for the first derivative with respect to x appears to be required to improve the accuracy. 

6.4 CORRELATION OF EXAMPLES FOR SWEPT WINGS OF VANISHING THICKNESS 

Calculations were performed at a Mach number of 0.9 and reduced frequency of 0.13 for sweep angles of 30 
deg, 40 deg, and 45 deg. Figures 41 through 43 show the real and imaginary parts of the jump in pressure coeffi- 
cient across the flat plate for the three sweep angles using a time step of 0.05. The greatest variation with sweep 
angle occurs at the plane of symmetry (figure 41). The solutions for a sweep angle of 30 deg are shown in figures 
44 through 46 for spanwise locations at fractions of semispan of 0., 0.51, and 0.93. The results after 150 and 300 
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iterations are compared with the integral equation solution of RHOIV. It appears that 300 iterations yield better 
results, and agreement with the RHOIV is seen to be excellent. 

The jump in pressure coefficient for 45 deg, the largest sweep angle for which the program converged, is 
compared with the RHOIV solution in figures 47 through 49. The agreement is seen to be as good as for the lower 
sweep angles. 

To test convergence at the higher frequencies, the 45-deg swept wing was calculated using k = 0.5. The 
solution converged as easily as for the lower frequency calculations. As seen from figures 50 through 52, the 
agreement with RHOIV is, however, not as good. As in the two dimensional solutions, the second-order dif- 
ference for the first derivative with respect to x appears to be required to improve the accuracy. 
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7.0 CONJUGATE GRADIENT TECHNIQUES 
FOR THE DIRECT SOLUTION 


7.1 DESCRIPTION OF BASIC CONJUGATE GRADIENT TECHNIQUES 

Although the direct solution as used in OPTRAN2 provides an efficient and practical procedure for analyz- 
ing two-dimensional configurations with mesh systems of the order of 4500, three-dimensional problems using 
mesh configurations of the order of 50,000 to 100,000 would appear to require alternate procedures to achieve 
comparable efficiencies. Review of the literature indicated that conjugate gradient techniques might provide a 
practical way to obtain a direct solution for larger systems. However, our problem is asymmetric and indefinite, 
and thus the classical conjugate gradient methods are not applicable. 

The classical conjugate gradient method for solving systems of linear equations of the form Ax = b whose 
coefficient matrices, A, are symmetric (Hermitian in the case of complex matrices) and positive definite was 
presented by Hestenes and Stiefel (ref 25) in 1952. This method converges to the true solution of the lineeu sys- 
tem in a finite number of iterations in the absence of round-off errors and hence may be thought of as a direct 
solution. However, the method was not widely used and little was heard about it until the mid 1960s. J.K. Reid 
(ref 26) noticed that the method is very efficient for large sparse symmetric positive definite systems that are 
well conditioned, with the asymptotic rate of convergence being inversely proportional to the square root of the 
condition number. In the last decade, many variants of the conjugate gradient method have been applied to large 
sparse problems with considerable success (See Hafez and Wong, ref 27). 

When the coefficient matrix, A, is asymmetric, the usual procedure is to multiply both sides of the equa- 
tions by the conjugate transpose, A*, to obtain a Hermitian positive definite linear system of the form 
A*Ax = A*b on which to apply the conjugate gradient method. However, the condition number of the matrix 
A* A is the square of that of A. That is, the asymptotic convergence rate of the conjugate gradient method applied 
to the above equation would be inversely proportional to the condition number of A instead of to the square root 
of the condition number of A. 

Various authors, such as Axelson (ref 28), Concus and Golub (ref 29), and Elman (ref 30), have generalized 
the conjugate gradient method to a larger class of matrices, so that the asymptotic convergence rate is inversely 
proportional to the square root of the condition number of A. Yet all of these methods still assume the coefficient 
matrix to have certain properties, such as A* -i- A being positive definite, and thus they are not always applica- 
ble to general nonsingular matrices. In appendix E, a variant of the conjugate gradient method is developed 
which is applicable to general nonsingular matrices and which has an asymptotic convergence rate inversely 
proportional to the square root of the condition number. This new algorithm together with examples of its 
application will be discussed in the following. 

However, for the new algorithm, when the coefficient matrix is even mildly ill-conditioned, the convergence 
rate of the conjugate gradient method is impractically slow. This difficulty, frequently encountered with the 
matrices for the transonic problem, can be remedied by preconditioning, which attempts to improve the condi- 
tion number of the matrix. Various authors, such as Lewis and Rehm (ref 31), Kershaw (ref 32), and Chandra 
(ref 33) have suggested different preconditioners for specific problems. The success of the method for a particular 
problem very much depends on the preconditioner chosen. For our problem we have tried preconditioners from 
the following sources: 

(1) A sparse capacitance matrix method 

(2) The ADI operator 

(3) The incomplete LU factorization procedure. 
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The preconditioner from the incomplete LU factorization has proved to be most efficient. We shall describe these 
three preconditioners in detail in section 7.2. 

7.2 PRINCIPLE OF PRECONDITIONING 

The principle of preconditioning is to recast the system of equations in the form 

(I + T)x = b 

where T is a matrix whose magnitude is small. One procedure for accomplishing this is to find a matrix N such 
that its inverse is easy to compute and N is very close to the coefficient matrix A. If such an N can be found, then 
the matrix 


R = A-N 

is of small magnitude. The matrix N*^R is also of small magnitude if A and N are well conditioned. We can write 
the equation Ax = bas(N + R)x = b. When we premultiply both sides of the latter equation by N'^, then we 
obtain 


(I - N“lR)x = N"lb 

Note that the singular values of (I + T) = (I-N'^R) (See Householder, ref 34) are defined as the square roots of the 
eigenvalues of 


(I+T)*(I+T) = I + T + T* + T*T. 

Since the magnitude of T is small, the singular values of the matrix (I + T) are close to 1. This is the ideal 
situation for applying conjugate gradient type methods. 

Aside from the case in which the matrix A is a so-called M-matrix, which is a matrix with negative diago- 
nal entries, non-negative off-diagonal entries with all its eigenvalues on the right half-plane, (see Varga, ref 35), 
there are no theoretical results in the literature concerning the choice of N. Almost all of the results in the 
literature concerning preconditioning are empirical and therefore are problem-dependent. The most popular 
and successful preconditioner is the incomplete factorization procedure. 

7.2.1 THE INCOMPLETE FACTORIZATION 

The concept of incomplete factorization may be used for preconditioning in the following manner: the L and 
U matrices resulting from the lower and upper triangular decompositions of the coefficient matrix (or some por- 
tions thereof) are modified to obtain matrices which are easy to invert and store. In the current algorithm, L and 
U are modified to be lower and upper tri-diagonal matrices. 

The incomplete factorization of the matrix A in comparison with the complete factorization of A is illus- 
trated in figures 53, 54, and 55. Figure 53 represents the sparsity of the original matrix A. Figure 54 shows the 
sparsity structure of the two matrices resulting from the complete LU factorization of matrix A. Figure 55 shows 
the sparsity structure of the incomplete LU factorization of the matrix A. The incomplete L and U are obtained 
from the regular Gaussian elimination by retaining only the nonzeroes in the sparsity structure of A. The 
details are expounded in appendix E. Our results are summarized as follows: 
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Problem 

Grid 

size 

No. of 

mesh points 

Formulation 

No. of 
iterations 

CRAY 
CPU sec 

(1) 

72x60 

2160 

ADKt = .3) 

263 

15.3 

(2) 

72x60 

2160 

direct 

346 

20.3 

(3) 

50x15x40 

15000 

direct 

564 

179. 


In comparison with other methods, the ADI described in section 6 takes 40 seconds for problem 1, and the 
out-of-core solver takes 4.5 seconds for problem 2, and 400 seconds for problem 3. 

In the two-dimensional problem, the ADI formulation (problem 1 of the test problems) seems to have better 
convergence properties than the direct solution formulation. This might well be true for the three-dimensional 
case also. 

As an iterative method, the conjugate gradient is not so “operator-sensitive/’ For instance, in problem 1 of 
the test problems, second order central differencing is used for the (p^t term. In the ADI method described in 
section 6, second order upwind differencing is used for the (p^t terms because the ADI diverges when central 
differencing is used. 

7.2.2 THE ADI OPERATOR AS A PRECONDITIONER 

The ADI method can be represented by the following matrix equation: 


<Pn+l 


(7.2.1) 


Assume equation (7 .2.1) converges, then = , <Pn = •Pn-i> equation 7.2.1 3 can be written as 

(I - G - P) = c. (7.2.2) 

We apply our conjugate gradient method (USYMLQ) to equation (7.2.1). Experimentally, we find that the rate of 
convergence of this method is similar to that of the ADI and hence more expensive than the ADI because for each 
iteration, we have to compute (I - G - P)x^ as well as (I - G - P)*x,^. 

We have only tried this on the 2D problems. Since this method is not comparable in efficiency with the ADI, 
we do not list the results. 

7.2.3 THE SPARSE CAPACITANCE MATRIX METHOD AS A PRECONDITIONER 

The method we discuss in this section is motivated by the two-step method discussed in Ehlers and Weath- 
erill (ref. 7). The flat plate equation without the boundary condition on the airfoil can be solved by Fast Poisson 
Solvers like FISHPACK (see Sweet and Swartztrauber, ref. 36). Since Ehlers and Weatherill (ref. 7) observed that 
the solution of the airfoil equation is close to the solution of the flat plate equation in the far field, we thought a 
matrix of the form: 


M = 



(7.2.3) 
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where S is the coefficient matrix of the difference equation of the airfoil equation for the near field mesh points 
(the mesh points inside the box which encloses the airfoil and all the supersonic points), and F is the coefficient 
matrix of the flat plate equation for the far field mesh points, would closely approximate the coefficient matrix of 
the airfoil equation (see fig. 56 below). In other words, if T is the coefficient matrix of the difference equation for 
the airfoil equation for the far field mesh points, that is, if 

- 

we assume (T - F) is of small magnitude. Therefore we use M in equation (7.2.3) as a preconditioner. Equations of 
the form M’^x = b are solved by a sparse capacitance matrix method. This preconditioner proves to be inefficient 
for the following reasons: 

1. The sparse capacitance matrix method requires an efficient complex sparse linear equation solver. The solver 
we were using, ME28A from Harwell, was written for IBM machines and is extremely inefficient for the 
CRAY. 

2. The two-step method discussed in reference 7 proves to be successful only for symmetric flow, thus the matrix 
M in equation (7.2.3) does not closely approximate the original coefficient matrix. Therefore convergence is 
slow. 

3. The box which encloses all the supersonic points can be almost as big as the full grid. 
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8.0 THE DIRECT SOLUTION FOR THREE-DIMENSIONAL WINGS 


The pilot program for the direct solution of two-dimensional configurations has been rewritten for applica- 
tion to three-dimensional rectangular configurations. This new pilot program, called OPTRAN3, is pro- 
grammed for the CRAY. Application of OPTRAN3 is discussed in section 8.1 and an improved algorithm for an 
out-of-core solver is discussed in section 8.2. 

8.1 APPLICATION OF OUT OF-CORE SOLVER 

One large and several small examples were run to test OPTRAN3. The large run, the only one to approxi- 
mate a practical configuration, was a 50 by 20 by 40 mesh system applied to an aspect ratio 3 rectangular wing of 
vanishing thickness. Since the flow for this configuration is symmetric top to bottom, this problem consisted of 
20,000 points. Pressure distributions from this calculation are compared with results from OPTRAD3 and 
RHOIV in figures 16 through 18. Correlation of all the results from the two finite difference procedures with 
reference results from the kernel function method is good and is about what we would expect from our experience 
with two-dimensional calculations. The overwhelming characteristic for the example is its running time which 
was some 1200 CPU seconds. Perhaps of even more significance is the fact that using the company cost 
algorithm, the CPU seconds only accounted for about one-fourth of the total cost of the run, while the remainder 
of the cost is mainly due to input/output operations. It is also noted that costs rise rapidly with increasing mesh 
points. For example, a problem with three-fourths the mesh points requires about one-third as many CPU 
seconds. 

The pilot program that was tested was developed on a version of the CRAY which has one million words of 
core storage. It was also developed during the installation period at Boeing of both the CRAY hardware and 
software. Now the operating system has been stabilized, and the original machine has been replaced by one with 
two million words of core storage. The larger memory allows us to bring larger blocks into memory. Larger 
blocks mean that the length of the vectorizable do-loops is longer, and execution time should be reduced signifi- 
cantly. It also means a reduction in the number of disk access operations. With the new algorithm, the number of 
words written on disk is also significantly reduced. When solutions for a number of right hand sides are required, 
the out-of-core solver may be competitive with iterative methods, such as the ADI and the preconditioned conju- 
gate gradient procedure. 

The total cost of the initial three-dimensional run using the direct solution program was large enough that 
it was decided to concentrate efforts on the ADI procedure of section 6. If the running time can be reduced signifi- 
cantly as described above, the advantage of being able to include a number of mode shapes (i.e., a number of right 
hand sides to the set of difference equations) with little increase in cost would make the direct solution again 
competitive with other procedures such as the ADI method. 

8.2 AN IMPROVED OUT-OF-CORE SOLVER PROCEDURE 

In the 2D program, the out-of-core direct solver, ETCSM, is a general-purpose banded solver. Because of the 
size of the 3D problem, we intend to design a direct solver specialized for our problem. The key idea in the new 
solver is “implicit factorization”, which means that the LU factorization of the diagonal blocks, which are prod- 
ucts of sparse matrices, is computed as needed. The “implicit factorization” requires the matrix to be block- 
tridiagonal. The following 3 by 3 block system serves to explain the difference between the algorithm used in 
ETCSM (the out-of-core solver of ref 18) and implicit factorization: 
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We are interested in solving Ax = b, where A is represented by the block system: 



”Aii 

Ai2 

0 


■bf 

A = 

A 21 

A 22 

CO 

b = 

b2 


0 

A 32 

As3_ 


CO 

1 


The matrix A is now decomposed into a product of block-lower and block-upper triangular matrices: 



-I 

0 

0" 


■Ull 

U12 

0 

A = LU = 

^21 

I 

0 

♦ 

0 

U22 

U23 


0 

L32 

I 


0 

0 

1 

CO 

CO 

ID 


where, by comparison with equation (8.2.1), we see that 
Uii = Ajj, Ui2 = Ai2 

L 21 = A2iA^j“^, U 22 = A 22 - ^21^12’ U 23 ^ A 23 
^2 "" A32U22"^ 

U33 "" A33 - L32U23. 


Thus Ax = bis to be solved in the form LUx = b.Lety = Ux.ThenLy = b can be solved by forward substitution: 

Yl = bi 

Y2 = b2 - L2iyi 

^3 " ^3 “ ^2^2- (8.2.3) 

The solution vector x, which satisfies the equation Ux = y, can now be found by backward substitution: 

^3 ■ U33"^3 

^2 = U 22 "^(y 2 ■ U 23 X 3 ) 

xj = Uifl(yi - U 12 X 2 ) (8.2.4) 
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In practice we do not need to compute the inverses of U33, U22, and U^, but can use their lower and upper tri- 
angular factors to obtain X3, Xg, and x^ by forward and backward substitution. Continued application of LU 
decomposition results in inverses being required only for diagonal matrices, a simplification that, in turn, 
requires storing only the original block matrices. Thus, for L21 = A2iVu-i, we store A21 and = Uy) instead of 
^ 21 - 

In the explicit factorization, as in the algorithm used in ETCSM, all the nonzero blocks of the factorization 
(except the identity matrices of L) are stored on disk, and the sequence of calculation is exactly as those described 
by equations ( 8 . 2 . 2 ) to ( 8 . 2 . 4 ). In the implicit factorization, we store only the diagonal blocks Uu, U22, and U33 in 
their factored form on disk in order to conserve disk space. The sequence of calculation for the factorization stage 
is the same as (8.2.2); the sequence of calculation for the forward substitution stage is different: in place of L21, we 
use its equivalent, A21A11-1; in place of L32, we use its equivalent, A32U22'^ (see equation (8.2.2)). 

Thus equation ( 8 . 2 . 3 ) becomes: 


yi 

= bi 


^2 

= b2~ A2iAii"lyi 


ys 

“ ^3 ■ -'^32-^22 ^y2 

(8.2.5) 


The sequence of calculation for the backward substitution is the same as described by equation ( 8 . 2 . 4 ), with the 

exception that we take advantage ofthe fact that U12 = Ai2,U23 = A23, and substitute Aj2 in place ofUi2,A23 in 
place of U23 in equation ( 8 . 2 . 4 ), which now becomes: 

^3 = U33~ly3 

^2 = U 22 \y2“-^23^3) 

^1 ~ ^11 ^(yj - A]^2^2) (8.2.6) 

In this approach, we do not need to store the off-diagonal blocks of the lower triangular matrix L on disk, and we 
can obtain the off-diagonal blocks of the upper triangular matrix U readily from the original matrix A. 

We have tried this approach on the two-dimensional problem; the cost is about one-third of that of ETCSM. 
For the 3 D problem, the matrix needs to be reordered with one level of nested dissection to reduce the bandwidth 
in order to use this approach more effectively. 
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9.0 CONCLUSIONS 


This investigation has centered on the development of procedures for calculating the unsteady transonic 
flow about harmonically oscillating three-dimensional wings. The work has included studies of the direct solu- 
tion method, a relaxation procedure using ADI techniques, and preconditioned conjugate gradient procedures. 

The most practical procedure for three-dimensional analyses requiring a limited number of mode shapes 
and using a computer with one million words of memory or less appears to be the ADI procedure. A pilot pro- 
gram, including a coordinate transformation for swept and tapered wings, has been developed and applied to 
several wings. Results appear generally reasonable and in good agreement with experiment and kernel function 
methods. Solution convergence characteristics appear similar to those of XTRAN3S, and results have been 
obtained for a swept, untapered wing with a leading edge sweep of 45 deg. However, it should be noted that the 
solution time is proportional to the number of modes. 

A pilot program for rectangular wings using the direct solution has been developed for the CRAY com- 
puter. As expected, this program has proved expensive to run for practical problems. However, a significant part 
of this expense is due to the cost algorithm that penalized input/output operations. This was a problem with the 
one-million- word CRAY on which the current program was developed. With CRAYs having larger core storage 
now available (we are currently working with a two-million- word core, and there is talk of an eight- to thirty-two- 
million- word core capability in the future), the direct solution remains a viable procedure, particularly for prob- 
lems involving a large number of modes, since solutions for additional modes are obtained at minimal extra cost. 

Finally, an algorithm for the conjugate gradient procedure as applied to asymmetric, indefinite coefficient 
matrices has been developed with a solution convergence rate proportional to the square root of the condition 
number. Several preconditioning procedures have been tested, and one based on incomplete LU decomposition 
has proved most efficient. Running times appear to be generally comparable to those for the ADI procedure. This 
procedure is still under development. 
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APPENDIX A 


DERIVATION OF THE ADI DIFFERENCE EQUATIONS FOR THE 

RECTANGULAR WING 

A.1 THE BASIC TIME-DEPENDENT DIFFERENCE EQUATIONS 

In the same manner as for the two dimensional ADI method, we consider the time-dependent linear equa- 
tion for the perturbation unsteady potential. Thus 

(*^tt ^ 2</>xt) h = («‘^x)x ^ ^ 

Following the procedure used by Borland, Rizzetta, and Yoshihara (ref. 17), we write the difference equation for 
the three sweeps as 

X sweep: 

26^ At) = ^ </>" + 

y sweep: 

25^ (/-^“)/(e At) = 6yy(</.^-<A”)/2 

z sweep: 

*"■!)/(. At^) . 2« At) = 

X 

The superscripts a and \ denote intermediate steps between the n and n + 1 iterations. The difference oper- 
ator is a backward difference, while is a central difference for u positive (elliptic) and backward dif- 

ference for u negative (hyperbolic). 

As in the two-dimensional version, we introduce harmonic motion and write for the nth and the next inter- 
mediate approximation 

, ^iigi o> n At^ ^ ^X^i a, (n+1) At_ ^ co (n+1) At 
Substituting these expressions into equations (A-2) through (A-4) yields 

( <p“- <^*^) /(e At) = 5jj(u6jj<p") /2 + + «yyV’" + ^^-5) 

26x(/-<^“)/(« At) = (A-6) 

- 2p^<p^ + /3iV’^)/(e At) + - <p^)/U At) = (A-7) 

where 

=e~ i At 


(A-2) 


(A-3) 


(A4) 
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To write the equation for coding, we use the same operators as in the direct two-dimensional method. For 
the difference operators, we write 

26^V>/(6 At) = jj^) 

6x(u5x<p)= 2cj - <^ijk)-2d. jj^) 

(Vljk" “ijk)>0’ 

2Ci.i ujjk (v’ijk - v’M jk)- 2di.i u-.j • - <p-^2 jk) 

(Vljk"“ijk)<0- 

6yy<^ = 2aj(<pj j.i k ■ ‘^ijk)" 2bj(v>ijk ■ ‘^i j+1 k) 

- <p^^y 2bj^(v,..j^ - <p.. (A-8) 


where C^i = 2/A t(xj - and the other coefficients are defined in reference 1. 


or 


A.2 THE X SWEEP DIFFERENCE EQUATIONS 

With the aid of equations (A-8), equation (A-5) for the x sweep becomes 

2^^<p“/(£ At) - 5 x(u6^<p“)/ 2 = At) + 5x('^^x‘^")/^ ^ ^ 

"3i(%k"^-ljk ) " ‘^i"i+ljk(^+ljk'^jk ) 

" ‘^i^jk(^jk ■ Vl jk ) = ^l[ "3i( %k ■ %-l jk ) 
jk(*^ i+1 jk *^ijk ) ijk i-1 jk ) 


(A-9) 


(A-10) 


Writing the equation in the form 


SUBl(I) * + SUB(I) * + DIAG(I) * , + SUPER(I) ♦ = RHS (I) (A-U) 

1-2 jk 1-1 jk ijk 1+1 jk 
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we see that 


SUBl(I) = 0.0 
SUB(I)= -C3i 

DIAG(I) = C3i + CiUj,jjj^ +d.u.j^ 
SUPER(I) = -CiUi,^.^ 

Note that DIAG<I) = - SUBfl) - SUPERfl). The right hand side term becomes 


where 


RHS(I) = j3j(RHSl + RHS2 + 2*RHS3 + 2*RHS4) + 2RHS4) 

RHS2 = CiUi, j ^ ) 

RHS3 = a. ^ ) - b. f <p”., - ) 

JV i]-lk i]k / ] \ ijk ij+lk/ 


When the point ijk is supersonic, then equation (A-10) becomes 

'3i (*’'ijk” ''i.ljk) " 'H“ijk(''ijk’ ‘'i-ljk ) * jk (*’° 1 jk " ‘'’”i-2ik ) 

['3i (■'“jk • ■'”.1 jk ) * 'U“iik (•'“jk ■ Vl ik ) - 'ii.l'H.l jk «.i jk ■ A-2 jk 

*^(-U.lk-%k)-^^iKk--"Hk) 

* 2“k k-1 ■ ■'"jk) ■ ^'’k (■'"jk ■ ■'"j ktl ) ] 

Comparison of equation (A-15) with equation (A-11) yields 


SUB(I) = -C3i + Ci.iUi.j^ + di.jU..j.j^ 
DIAG(I) = C3i-Ci.iUi.j^ 

SUPER (I) = 0.0 


(A-12) 


(A-13) 


(A-14) 


) 


(A-15) 


(A-16) 
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Note that SUB(I) = - SUBKI) - DIAGKI). The terms RHSl and RHS2 become 


RHS2 = (v” j - ^1.2 j ^ ) (A.17) 

RHS3 and RHS4 remain unchanged. 

A.3 BOUNDARY CONDITIONS FOR THE X SWEEP 

We apply the outgoing wave type boundary conditions of Engquist and Majda (ref, 11) on the upstream and 
downstream boundaries. Following the two-dimensional method, we obtain 


■'“jk ' °kl"2jk * (*'2jk ■ °kl‘'ljk ) 

‘m.vik " ik * imax' Jk) 


max-' 


where 


^kl ^ M(x 2 - xi)/[(l- M) At] 
Cki - (1 - Cki)/(1 - Cki) 


(A-18) 


(A-19) 


Ck3 = (l-Ck3)/(l^Ck3) 


(A-20) 


The equations (A-12) are modified for i = 2 and i = i^ax'l 

DIAG(2) + DIAG(2) + SUB(2)=kCj^l 

RHS(2) = RHS(2) - SUB(2)*/?i(<^^.j^ - 

SUB(2) = 0.0 (A-21) 


DIAG(IMAXl) = DIAG(IMAXl) + SUPER (IMAXl)*Cj^3 
RHS(IMAXl) = RHS(IMAXl)-SUPER(IMAXl)*i3, 

\ ^max 

SUPER(IMAXl) = 0.0 


-1 j k imaxi^) 


(A-22) 
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For k = just below the wing plane, and k = k^ + 1, we apply the wing boundary conditions when 
i^ < i < i^, for j < jg, the index of the first spanwise variable beyond the wing tip. The conditions are in the form 

(9v’/9z)z - - 0 = F^(x,y) 

(9<P/9z)^_^Q = F^(x,y) 

For k = k^, we have 


where h = Zj^ ^ Similarly, for k = k^„ + 1, 

^km+l(’’ij V km+l) ' " ^^k„i+l ^ ij 

We see that for k = k,^, the right hand side term is modified by 

RHS(I) = RHS(I) + 2^,b, [hpL +^n - u 

1 ^mL ij ijkjn iJ k^+l 


(A-23) 


(A-24) 


(A-25) 


Similarly, for k = k^+lwehave 


RHS(I) = RHS(I) - 2a, 


-uuU . n n 

nb 

. iJ ijkm ij kjn+1 


(A-26) 


For i > ii and j < jg, we must satisfy the continuity of the normal derivative across the wake sheet and 
continuity of pressure. Detailed discussion of this is given in reference 1. The jump in the potential at 
X = Xj^ + 1 is chosen to satisfy the Kutta condition that the jump in pressure at the trailing edge be zero. Con- 
tinuity of pressure is assumed by setting 

= Vl) 

where Acp^j denotes the jump in potential at x = Xj. The continuity of normal derivative is assured by adding to 
the (p 2 z term the quantity 


- '>k„ ^'^ij 

for k = kj^ and 

^km+l ^‘^ij 

fork = k^ + 1. 
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A.4 THE DIFFERENCE EQUATIONS FOR THE Y SWEEP 


The difference equation for the y sweep, equation (A-6), becomes 


or 


- 2E^<p^/(,e At) + 5yy//2 = i3j6yy<p“/2 - At) 


M k - »yk) ■ jn k) ■ '3i(%k ■ M jk) 

= j.i k' •'"ik) ■ ’’)(■' «k ■ '’1 jti k)] 

- '3i(*“jk ■ <1 jk) 


Writing this equation as 


SUB(J)* , + DIAG(J)* + SUPER(J) * ip\ . , , = RHS(J) 

1 ]-l k iJK 1 j+i k 


and comparing equations (A-28) and (A-29), we obtain 


(A-27) 


(A-28) 


(A-29) 


SUB(J) = aj 

DIAG(J) = - aj - bj - cgj 
SUPER(J) - bj 

RHS(J) = j ^ J 

■'=3i«jk‘ Vljk^Vljk) 

Note that (p^.i j ^ is known from the previous sweep, since we always move in increasing i. Since the (p^ and cp“ 
values are not saved, the quantity 

a 

^ i-1 jk 


has been written over with (Pu jk by the previous step in the y sweep. Hence, after the i - 1 sweep we must save 
<Phjkby 


PHIOLD(J) =<p“ ^ (A-31) 

before storing the results of the i - 1 step. In each step of the y sweep the value of k is fixed, while i varies from 
i = 2 to i„ax ■ 1- Then k is increased by 1 and i is varied from 2 to i„^ax • 1 antil the entire grid is swept. 
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A.5 BOUNDARY CONDITIONS FOR THE Y SWEEP 


The difference equations for the y sweep require boundary conditions on the upstream x boundary. Hence, 
the value of i = 2 is a special case. Then the term “CgjCPj.i is missing from the right hand side RHS(J) term, and 
we replace it with the expression 

‘"3i^kl[‘^2jk' ‘^2jk] 


The diagonal term is modified by 


DIAG(J) = DIAG(J) + CggCj^^ 

The line y = 0(j = 2) is a plane ofsymmetry for the solution since it is the location ofthe root chord. Thus 6 
(p/8y = 0 at y = 0 yields 


*^ilk ‘^i3k 

Application of the boundary condition to equation (A-29) yields 


SUPER(2) = SUPER(2) + SUB(2) 
SUB(2) = 0.0 


(A-32) 


(A-33) 


We apply the Majda and Engquist nonreflecting boundary conditions at the outer spanwise boundary, 
y = (yj max + yjmax - 1^2. This takes the form 


i imaxk ^^2*^ i k i jmax l k ^k2‘^ i k) 


(A-34) 


where 


^k2"(^"^k2)/(^‘^^k2) 

Ck 2 = M/'yj -y: VIT/ (l - M2)At 

\ Jmax Jmax 

K = (l - M^)/(M^e) 


(A-35) 


This requires the following modiflcations 


DIAG(JMAXl) = DIAG(JMAXl) + Cj^2*SUPER(JMAXl) 

RHS(JMAXl) =RHS(JMAX1) - SUPER( JMAXl) * /3, ^ v’". , , ,\ 

Umax-1 k ijmaxkj 

SUPER(JMAXl) = 0.0 (A-36) 
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A.6 THE DIFFERENCE EQUATIONS FOR THE Z SWEEP 

The difference equation for the z sweep in equation (A-7) may be written as 

V2 - At^) - V(e At) 


= - {2<p^ - )3iv>''-^)/(€At2)] - 2^//(e At) 


With the aid of equation (A-8) we obtain 
n+1 n+1 


- -Sk ) - - C:) - "i<k- 

° « k-1 ■ ‘'“jk) ■ '’k(%k ■ « ktl) 


Writing equation (A-38) in the form 


SUB(K) * , + DIAG(K) * + SUPER(K) * cpVi\^ , = RHS(K) 

ij k-1 ijk ij k+1 


and comparing with equation (A-38) yields the following relations for the coefficients: 
SUB(K) = aj^ 

DIAG(K) = - - bj^ - - Cgj 

SUPER(K) = bj^ 


RHS(K) = ^ ^ J J 

■ ■ '*l''!jk) ] ■ '3i(4k - U jk* ■'“l jk ) 


(A-37) 


(A-38) 


(A-39) 


(A40) 


where Ej = 1 /sAt^. Note that the right hand side term, RHS(K), contains the known value of the n - 1 - 1 approx- 
imation, Since the n + 1 approximation replaces the X approximation, (p^j j ^ must be saved in the pre- 

vious step. Thus we define 


PHIOLD(K) = 

1-1 jk 


(A41) 
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A.7 BOUNDARY CONDITIONS FOR THE Z SWEEP 


For the z sweep we must apply the boundary conditions on the wing, on the wake, and on the mesh bound- 
aries except for the side boundaries, that is, the x-z plane boundaries. On the upstream boimdary, we use 




n +1 

Ijk 



(A42) 


No boundary conditions are actually needed in the difference equations on the downstream boundary, but for the 
next step we need the values of (p^ which satisfy the downstream boundary conditions. Thus we set 


where 




n +1 _ 7 ^ n +1 


^maxil^ 


jk ^ ^l( imax-1 ik imax jk) 


^i.q = M/x- -X- ,\/(M+l)At 

\ ‘max ‘max'‘/ 


(A43) 


^k3" “ ^k3)/(^ ■^^ks) 


(A-44) 


On the upper boundary, we have for the nonreflecting boundary conditions 

n +1 n +1 o / _n _ 7 ^ .n ^ 


_T; n +1 X o / n 7 ; n 

^max ^max“^ ' ^max* 


k .1 ^kS^iik 


(A45) 


max' 


where 




(^k 


■ ^k 

max ^max 


■1) 


vie /[(l- M^)At] 


Ck5 = (l-Ck5)/(l^Ck5) 


(A46) 


The lower boundary nonreflecting condition is 

n+ 1 - 7 ^ n+ 1 , n / n 7 ^ n \ 

"’ijl ■^k4‘^ij2 "^l(‘"ij2‘^k4<^ijJ 

where 

Ck4 = M(z 2 - Zj) y/K/{\ - M^)aI 


(A47) 


Ck4 = (l-Ck4)/(l^Ck4) 


(A48) 


On the wing, for io < i < ii where io and i^ may depend on j, we have for z = (zj^^ + +i)/2 = 0, 

(d(p/dz)^ = F^. 

ij 

(9.^/dz)'' = 

ij 
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For k = , the term ^ i) becomes 


b. +i'\=-bu hai 2 _=-b. HFV; 

^m\ ijkm y dz kjjj ij 


(A49) 


Fork = + 1, the term akm(<P"jk.i - (pPj,J- ^ ) becomes 

^n +1 n +1 

ijkm ij •'m 




u 


kin+l) " ^ ^km+1 ^ ij 


(A-50) 


These same boundary conditions must be applied to the corresponding terms for cp" as well. On the wake for i > ij, 
to satisfy the continuity of the normal derivative across the wake, we must add for k = k„ the terms 


5], (<P-- 

•^mV ij 


n +1 n +1 n +1 . 

■ 1 . ~ 


ijkm ijkm+1 L '"U 
where Acp = (py - (p^. Similarly, for k = k,^ + 1 we have 


n+1 ) = bi / \ + fi A ”■'■1 

U ' ^m\ ijkm ^ijkm+l) ^km^^^ij 


(A-51) 


^ / n +1 n +1 \ ^ A n +1 

^■‘■1 (^ij kjn ^i} kjjj+l) 


(A-52) 


Equations (A-51) and (A-52) are applied to cp" as well. 

For all k, the boundary conditions on the upstream boundary require the following modifications of the coeffi- 
cients. When i = 2, 

DIAG(K) = DIAG(K) + C 32 

RHS(K) = RHS(K) - C 32 ~ <^ijk ^^ 1 ) (A-53) 

On the lower boundary, for k = 2 

DIAG(2) = DIAG(2) + Cj^4*SUB(2) 

RHS(2) = RHS(2) - SUB(2) *%( v’”j 2 - 0^ 

SUB(2) = 0.0 (A-54) 

On the upper boundary, for k = k„,ox ■ 1 = KMAXl 

DIAG(KMAXl) = DIAG(KMAXl) + SUPER(KMAXl)*Cj ^5 

RHS(KMAXl) = RHS(KMAXl) - SUPER(KMAXl)*)3i f , - C, , ) 

okmax -1 0 kmax' 

SUPER(KMAXl) = 0.0 (A- 55 ) 
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For the wing boundary conditions at k = k,„, we have 


DIAG(KM) = DIAG(KM) + b. 

SUPER(K) = 0.0 

RHS(KM) = RHS(KM) - b. hF^ + )8,rb. (/\., - , ,Wb. hF^.l (A-56) 

Km ij IL ^mV ijkm U km+l^ yj 


For the wing boundary conditions at k = + 1, we have 


DIAG(KMP) = DIAG(KMP) + a, . , 

Km ^ 

SUB(KMP) = 0.0 


RHS(KMP) = RHS(KMP) - /ii[h pU t 




n 


ijkm+1 




For the wake boundary conditions, at k = k^^ we have 


RHS(KM) = RHS(KM) + 

RHS(KMP) = RHS(KMP) - a, - p, Av>".') (A-58) 

Km'^J^V ij ^ ij/ 


A.8 PROVISION FOR THE SYMMETRY OF THE STEADY FLOW 

When the steady flow is symmetric and 

pU = pL 
ij ij 

then the unsteady potential is antisymmetric, and we have 

‘^ijkm+1 ^ ■ ‘^ijkm 
*^ijkm''’2 *^ijkm‘i 


(A-59) 


(A-60) 


These may be employed as boundary conditions, and half the matrix can be solved with a considerable saving in 
computing cost. 

In place of the boundary conditions for K = KMAXl -1), we apply boundary conditions for k = k^. 
Thus the asymmetry conditions yield for k = k^^ and for i > \q and i < ij, 

DIAG(KM) = DIAG(KM) -b, 

SUPER(KM) = 0.0 (A-61) 
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A.9 EXACT EQUATION SOLVED BY THE ADI METHOD 

The ADI method of relaxation leads to a complex potential in much the same way as the block relaxation 
procedure that is used to solve the harmonic oscillation equation. The time step At is irrelevant to the solution, 
but its use leads to truncation errors which are of the order of At. To find the differential equation solved by the 
ADI method, we add equations (A-5), (A-6), and (A-7). We obtain 

At2) - At) 

+ (A-62) 


Let(p"'*'i = (p" = (p, then the differential equation can be recognized as 


(^‘^x)x ^ *^yy ^ ^zz 


1 + %) eAt 


- hf 


<p = 0 


(A-63) 


Letting At— 0 yields the classic unsteady equation for the perturbation potential. 
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APPENDIX B 


DERIVATION OF THE ADI DIFFERENCE EQUATIONS FOR 
SWEPT AND TAPERED WINGS 

B.1 BASIC COORDINATE TRANSFORMATION 

In analyzing swept wings, it is convenient to use a coordinate transformation in which the leading edge 
and trailing edge are defined hy single values of the streamwise coordinate. The partial differential equation for 
time-dependent transonic small pertxu-bation flow is given by 

(2<^xt ^ ‘^tt)/‘ = ("‘^x)x ^ ^yy ^ 

This equation is the classical form of the transonic small perturbation equation and does not include the higher 
order sv^ept \ving terms. (To obtain the basic steady state potential, cpo these higher order terms must be included 
to locate correctly the shock wave.) 

A coordinate transformation, which meets the foregoing criteria, is shown in figure B-1 and is described by the 
following equation, 


f = [x-XLE(y)]/S(y)- 1 
v = y 

I = z (B-2) 


where S(y) is the semichord of the wing cross section at the spanwise location y. Then 

</>2 = (B-3) 

The differential equation (B-1) becomes 

(20^/8 + (u + (B4) 

Multiplying by S yields 

(2<A^t + S%)/e = (u ^y ^ ^ ^ ^ 

The second and third terms from the last can be combined to obtain a conservation form of the differential equa- 
tion. Thus 
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For convenience, we define the following quantities 




(B-7) 


Y = SVy 

(B-8) 

Then the differential equation becomes 





(B-9) 


Following Borland, Rizzetta, and Yoshihara (ref. 17) we construct an ADI procedure for the solution of this differ- 
ential equation. We obtain 

For the ^ sweep: 

) = 6f[(u/S)((A" + 4>^)/ 2 ] 

+ 8 (S8 <t>^) + + 6 Y” (B-10) 

»7\ ^7 / K f 

For the Tj sweep: 

[(2/(e At)]^f {<l>^ - <#.") = ]/2 (B-U) 

For the sweep: 

[(2/(e At2)] - 20 " + + [2/(e At)] 8^{<f>^^^ - 0^ = - «/«") (B-12) 

As in the rectangular wing, we represent the nth approximation and the intermediate approximations between 
the nth and n + 1st in the form 

.n ^ i n o> At n , A _ -i(n+l) w At A 

,a _ J(n+1) o) At a 
<p - e ^ ' ip 

Substituting into equations (B-10), (B-H), and (B-12), we obtain the following equations: 

For the ^ sweep: 

[2/(e At)]'6^(v“ - = 6^[(u/S)6^(<p" + /3i<p”)]/2 

+ + /3j6^Y" (B-13) 
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where X” = + 6^(p*^) and = S6^((p*^^^). Writing the equation with (p« on the left hand side yields 

[2/(6 At)]^^<^« - 5^[(u/S) 5^^«]/2 = /3^ |[2/(6 At)] 

+ 6f[(u/S) 6^v’»^]/2 + + S6^^<pn + (B-14) 

For the q sweep: 

[2/(6 J /2 

Writing the equation with cp^ on the left hand side yields 

[2/(6 At)]‘5^vj‘'^ - 8^(^S8^,p>^y2 = [2/(e AtllV^.p® - p^(^S8^<p^y2 (B-15) 

For the C sweep: 

[S/(6 At)] (<pn+l- 2)3j<^n + ^2^n-l^ + [2/(t At)] 6^(<^n+l _ ^k) = Sg^^^^n+l _ ^^^n^ 

writing the equation with (p*^ ^ on the left hand side yields, after dividing by S, 

5^^^n+l/2 - ^n+1/ ^t2) _ [2/{e At S)] 

= p^^8^^<p^/2 - [l/(e At2)]^2<^n - j - [2S/(6 At)]^<^^ (B-16) 

In the ^ sweep, when we replace u by u = u/S, then the terms 6^ (u5^) are the same form as the unswept 
coordinate system. Similarly, 

S,(S6,.,'')/2 = ajSj. j ^ - „■!. J - bjS., ./, ( J (B-17) 

In the ^ and r\ sweeps, replacing 


®jSj. 1/2 = aj 


^jSj+ 1/2 = bj 


the term 5^(S8^(p") becomes 


^ ^j(‘^"i j-1 k ■ ‘^'ijk) ■ 2 bj( ■ ‘^"i j+1 k) 


(B-18) 


(B-19) 
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The difference equations for the ^ sweep become 


'=3i(‘'“jk ■ ■'“.1 jk) ■ jk ■ %k) * “jk ■ jk) 

° ^l['3i(» yk - 1 jk) * )(■'"♦! jk - ‘'"jk) 

■ ‘*i®y(‘'"jk ■ 1 jk) * ^“j(‘'" j.i k ■ •'"jk) 

- I,) * 2a^Sj(«'!. ^ J 

-^»k(%k'''"ikH)*¥"*V"] 


(B-20) 


where Cg^ = 2/[e At(Xj - Xj^)]. The difference equations for the ^ sweep of the swept wing differ in form from the 
rectangular version only in the factor Sj multiplying and and the addition of the cross terms, 


6^X^.5^Y^ 


We require the difference form of the derivatives, 6^X*^ + for the right hand side of the equations. 

With GK^,ti) = we have 




(B-21) 


The first term is the same form as 5^(u5^(p*^), with G^ replacing u, and can be easily written down. The Y'' term is 


= 6^(SG6^<pn^ 

From reference 1, the second order difference for 6^ is seen to be given by 

V lyk " jk ■ ^ - “^i-l jk) 

where and d^^ are defined on page 40 of reference 1. Similarly, we write 

V " j+1 k ■ ‘^ijk) ^ ^lj(‘^ijk ■ <^1 j-1 k) 

where the subscript j denotes that c^j is defined in the same way as qj but with the variables . 
Combining the two difference equations, we obtain for the cross derivative term 

-k - (c jj - dj-)GjjHjjk - d^jGj.j jHj^j jk 


(B-22) 


(B-23) 


(B-24) 


(B-25) 


where = Cy (cPy + u, - (pyk) + dy (CPyk - (Pij.lk) 

Similar results may be written down for 8^[SG5^cp]. No special application of boundary conditions is required 
since the boundary values of (p are calculated from the boundary conditions at each sweep. 
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Similarly, we obtain for the t] sweep 


■ '3i(‘'U ■ ■'m jk) ' ■ '8i(»yk ■ ”1-1 jk) 

*^i[^(-iMk-"jk)-^'i(-«k-"iHk)] 

This is seen to have the same form as the rectangular wing version. 


Finally, the sweep becomes 

n+1 n+l\ 1 / n+1 n 


„ / n+1 n+l\ 1 / n+1 n \ i? n+1 

ky^^ijk-l ^ijk / k(^ijk ^ijk+l) l^ijk 

- ¥3i( V - 

= ^l[“k(''"i k-1 ■ ■'ijk) - *>k(’' yk - % k*l) 

-®i(%k-vs;^)]-¥3i(-u-<iik) 

This has the same form as the rectangular version except Cgj is replaced by SjCaj. 


(B-26) 


(B-27) 


B.2 WING AND WAKE BOUNDARY CONDITIONS 

The wing boundary conditions are of the general form 

</’2 = f (x) + i w f(x) 

where z = 6flx)e'‘”‘ is the motion of each cross section. In terms of the variable this becomes 


Since 


then 


<^2 ~ f (f)fx i 


f=(x-XLE)/S- 1 


V>2 = f'(f)/S + io,f(f) 


(B-28) 


This differs from the rectangular version in the term Sj dividing the slope f(^). For the boundary conditions on 
the wake, we have 


or 


+ i CO A(p = 0 


+ icoSA<p-0 


(B-29) 
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In the wake boundary conditions, the reduced frequency is replaced by the product of reduced frequency and Sj. 

B.3 MESH BOUNDARY CONDITIONS 
On the upstream boundary, the condition for outgoing plane waves is 

M0^/(1-M), = O 


In swept wing coordinates this becomes 




(B-30) 


This is the same form as the rectangular wing version except for the factor of the semichord on the quantity 
M/(l-M). Hence the parameter now depends upon j. We then have for the boundary conditions 


^Ijk klj‘^2jk ^l(‘^2jk 


n+1 


where 


C|; jj = ({2 - {p MSj/(l-M)At 
C,y = (l-C^i.)/(I.C,y) 
Similarly, for the downstream boundary 

+ MS0^ / (1+M) = 0 


This leads to 


where 


n+1 ^ p n+1 ■ „ / n p, n \ 

imaxik k3j jk Wx 

Ck3j = (‘-Ck3j)/(>*C^3j) 

B.4 CONDITIONS OF SYMMETRY AT THE ROOT CHORD PLANE 

At the root chord plane y = t| = 0, we apply the condition of symmetry given by 

ip =0 

In the swept wing coordinate system this becomes 


(B-31) 


(B-32) 


(B-33) 


(B-34) 
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WeletG{^,Ti) - and let j = 2 be the plane of symmetry y = 0. In difference form, the conditions of symmetry at 

the point Xj become 


°i2(*’^2j ■ 4l2j)/(fi ■ fi-l) * (•''iSk ■ ” ak)/(’3 ' "l) = « 

Here we use backward difference in ^ since we sweep in the direction of increasing i. Since r\i = -TI 3 , 

*^4i ^ ~ ^i-l) ^ ^’’3^i2/(^i " ^i-l) 

The boundary condition then becomes 


411*^ i2k ‘^i-2k) *^i3k ‘^ilk"° 

from which 

X-X „/X X \ 

For the outer boundary condition, we have 

d>y +V^ Md>^/ ( 1-M2) = 0 
In swept wing coordinate, this becomes ' 




We write equation (B-38) in implicit difference form, for j = - 1 . We obtain 




1 j+l k ijk 


»7j+l - 


J 


. n 

(l-M2)At^ 


n +1 _ 

j+ V 2 k ^j+ V 2 k 


)=0 


Let Cgj G, jmax-l(’’jmax " ’’imax'l) ^(^1 ~ ^i-l) 
and C ^2 = Mv1C(,j„^- ,j„^,,,)/[(l-M2)At] 

then the boundary conditions become 


^ / n +1 ^ n n +1 « n \ . n +1 . « n 

-2i(»ijk '>1'0 ijk n.l jk ■ H jk) •'i )>1 k i 


j+lk 


_ *^■'■1- « . f~, / n+1 , n+1 on o n \ _ ^ 

■'ijk ijk ‘=k2(''i j,i k * ■'ijk ■ '’l'' i j*l k ■ ^1'' ijk) • ® 


n+1 


n+1 


(B-35) 

we define 
(B-36) 

(B-37) 

(B-38) 

(B-39) 


(B40) 
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Rearranging the terms leads to 


n+1 _ /p, _ :r X n+1 _ n+1 

i Jmax k " ( k2 21)^^! k + ^2i%i k 


[^2i( i j- „x-l k i imax-2 k) ^k2*^ ij 


' max 


^ ^ imax 


o] 


(B41) 


where C2i = C2i/(1 + €^2) and 0^2 = (1-Ck2)/(1 + 0^2). 


The boundary conditions on the wing and the wake are treated in the same way as for the rectangular 
wing, with the additional factor of the semichord as described in the preceding section. The boundary 
conditions on the outer ^ mesh boundaries are unchanged from the rectangular wing since the ^ variable is 
unaffected by the swept wing transformation. 

B.5 DERIVATION OF THE SWEPT WING TRANSFORMATION 

To illustrate the method of deriving the swept wing coordinate transformation, we consider the simple swept 
tapered wingwhose leading and trailing edges are straightlines.Let0s„ be thesweptangleoftheleadingedge and 
R be the taper ratio. Then the x coordinate of the leading edge is given by 

Xle = - 1 +y tan (B42) 

For the trailing edge, we write 


Xte = 1 + by (B43) 

and determine b so that at y = yj (that is, at the tip), the chord of the wing is equal to 2R. We then obtain 

b = 2(R-l)/yt + tan dgxv = 2C2 + Cj (B44) 


where 


C2 = (R-l)/yt and = tan 6^^. (B45) 

Because of the slope singularity in the pressure at the wing tip, the tip is placed halfway between grid 
points in the spanwise direction. The grid lines which coincide with the leading and trailing edges of the wing 
are extended beyond the last spanwise grid point on the wing planform, j = js, using a quadratic to preserve 
continuity of slope. At a point midway between the wing tip and the intersection of the linear extensions of the 
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edges, y„,, the slopes of the edge coordinate lines are both set equal to the average of the edge slopes at the tip. 
The edge coordinate lines are extended beyond the point y„, in linear fashion. The point of intersection of the edge 
extensions is found by setting Xle = in equations (B42) and (B43). 

y = -l/2C2 (B46) 


Let yjn, be the grid value of y for which 


^jm+1 

Thus for j < jm, we define 

Xle = -1 + Ciy + bi(y - yjs)2 

Xte = 1 + (Cl + 2C2)y + b2(y - yjs)^ 

The coefficients and bg are determined so that at y = 

X'le = Xte = [Cl + (Cl + 2C2) ]/2 = Ci + C2 

This leads to two equations for bj and bg which yield 

^2 “ C 2 /2 (yjni ^js) ” “^1 

Beyond y = the trailing edge and leading edge are given by 

Xle = + Ciyjnj + bi(yjj„ - yjg)^ + (Ci + C2)(y - yjm) 

Xte = 1 + (Cl + 2C2)yjm + b2(yjm ■ yjs)^ + (Cl + C 2 )(y - yjm) 
In summary, we have, for 0 < y < yjg, 

Xle = “1 + Ciy 

Xte " 1 (Cl + 2C2)y 

S(y) = (Xte ~ Xle)/2 = 1 + C 2y 

foryj3 :sy ^ yj^. 

Xle = -1 + Ciy - C2(y - yjs)“/2(yjm - yjs) 

Xte = 1 + (Cl + 2C2)y + Cg(y - yjg)^/2(yjm - yjs) 


(B47) 


(B48) 


(B49) 

(B-50) 

(B-51) 


(B-52) 


fB-53) 



and for y > yj^, 

Xle = -1 + Ciyjm - C2(yjm - yjs)/2 + (Ci + C2)(y - yjm) (B-55) 

Xte = 1 + (Cl + 2C2)yjm + C2(yjm ~ yjs)/2 + (Ci + C2)(y - yjm) ®-56) 

S(y) = 1 + C2yjm + C2(yjm " yjs)/2 (B-57) 

The transformation is given by 

iiv) = [x - Xle(^)]/S(^) - 1 (B-58) 

The mapping function required at every grid point is 

G(^,r?) = = - S' (r?)(f + 1) - X' le/S(^) (B.59) 


where X'^e S' are found by differentiating the expressions in equations (B49) through (B-57)). 


B.6 ELIMINATION OF SLOPE DISCONTINUITY AT PLANE OF SYMMETRY 

Some discontinuity affects at the line of symmetry can be alleviated by making the lead and trailing edges 
bend normal to the plane of symmetry. For j > 4, the coordinates of the leading and trailing edges are described 
in the foregoing section. The portion of the leading edge for j = 2 to j = 4 is fitted to a quadratic of the form 

XLE=XLE(0) + ay2 (B-60) 

The quantity a is determined to match the slope at j = 4. Hence, 

a = Ci/ 2 y 4 


Thus at j = 2 or 3 we have 

Xle = Xle(O) + Ciy2/2y4 

We determine Xl^CO) so that at y = Xle matches the known value, Xle 4 - Thus 

^le(O) + Ciy4/2 = Xle4 


and we finally obtain 

Xle = Xle 4 " Ci(y4 - y2/y4)/2 

Similarly, 

Xte " Xte 4 ~ (Cl + 2C2) ^y4 ~ y^/y4^/2 


(B-61) 


(B-62) 
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We see that 


S(y) = S(y 4 ) - C 2 (y 4 - y^/y4)/2 (B-63) 

Equations (B-61) through (B-63) replace equations (B-49) through (B-51) for j = 2 and 3. For larger sweep angles, 
the region should contain more than j = 2toj = 4 to adequately represent the geometry in the differencing. 
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APPENDIX C 

REVISED DIFFERENCE OPERATORS 

Cl THE OPERATOR FOR A CHANGE IN TIME STEP 

The- differencing in the equation for the x and y sweeps remains unchanged when the time step At is 
changed. In the z sweep, the only term which is affected is the second derivative with respect to time. Let the 
time step between the n-1 and the nth approximation be Ato and between the n and n + 1st approximation be At^, 
We express the time derivative in difference form as 

+ h(f>^ + c<f>^ ^ = 0 ^^ 

and expand the equation about the nth approximation. We then determine a, b, and c so that the left hand side 
gives at the nth approximation. Thus we obtain 

2 2 
a[<^ + + - ] + b</. + c[</, - "■**•] = 


Equating coefficients on the right and left sides of the equation yields 

a + b + c = 0 


aAt^ - cAtQ = 0 


aAt? + cAt? = 2 


Solving the last two equations simultaneously yields 

a = 2/|^Ati^Ato+ Ati)] 

c = 2/[Ato(Ato + Ati)] 

from which we obtain 


b 


(AtQ + Atij 


Introducing the variable R = At^/AtoandEo = (1 + R)/2R, E^ = HeEoAti^]; then 

6tt0/« = - 2EoRd>^ + Rd>"'^] 
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C.2 SECOND ORDER BACKWARD DIFFERENCE FOR 

For the first derivative with respect to x at the point X;, we write 

a(<^ij - = v’x 

Expanding about the point Xj, we get 

Ax^ 

- ip^ *** j 


- bj^<p- +_3 v?xx <p + (Ax ^ + Ax 2 )(Px - (^^1 ^ ^^2) v»xx '*'***]■ 

where Ax^ = Xj -x^,^ and Axg = x^.^ - Xj. 2 . Setting the coefficient of on the left hand side equal to unity and the 
coefficient of cp^x equal to zero yields 

a Ax^ - b Ax 2 = 1 

a Ax^ + bj^Ax^ - (^^1 ® 

We solve the two equations simultaneously for a and b and obtain 

a = ^2 Ax^ + Ax 2 ^/[^Axj^^Ax^ + AX 2 ) j 

b = Ax^/ Ax 

Writing 

^3i(^ij ” j) “ j “ ^i-2 j) - ^x 



yields finally 

Cgi = [2(xj - Xj.i) + Xj.i - Xi.2]/[(Xi - Xj.2)(Xi - x^.j)] 
^3i = (^i - *i-l)/[(*i - *i-2)(^i-l - ^i-2)] 
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APPENDIX D 

SHOCK BOUNDARY CONDITIONS FOR THE ADI METHOD 

Shock boundary conditions were derived by Hafez, Riszk, and Murman (ref. 37), by Williams (ref. 38), and 
by Seebass, Yu, and Fung (ref 39). The results were reported in reference 6 and extended to harmonic motion for 
application to the direct solution of the difference equations. The condition that the potential be continuous 
across the shock is given by 


[c^] = 0 

where [ ] denotes the jump in the quantity across the shock. The condition of continuity of mass across the shock 
from equation (A-11) of reference 7 is 

^-^ + <k-(y+l)> = 0 (D-1) 

where < > denotes mean value of the quantity at the shock. We now expand the shock conditions about the 
steady location Xq, with the shock position at time t being given by Xq + aX^. Then continuity of the potential 
across the shock becomes 

M = [cpq + + -] 

<^0] = 0 (D-2) 

= 0 (D-3) 


Similarly, expanding equation (D-1) yields 
2a ax, 

- ^^ + < k - (6 + 1) + a<pj^ + aXj(cpQ^^ + ' * ]> = 0 

Equating the coefficients of each power of a equal to zero yields 

C^k ~ ® 

2 i9X 

Tar •'lx ''Voxx> = “ 

Eliminating Xj by equation (D-3) yields 

! ^ i)[‘^ox]< ■'lx > - (ri) < •foxx > ['’1 1 = » 

Since K - (y + l)(Pox = u, we have 


(D4) 


(D-5) 

(D-6) 


(D-7) 


(D-8) 
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In implicit difference form, this becomes 


♦ <“x> [*’l]/2 - 1“1 < ^lx>/2 

In terms of the variables on the grid, we write 

Co3 = [u] = Ui^l.-Uij 

Cqj = 2/(€ At) 

Cq 2 = < «x> 


From the equation following (A-22) of reference 7, we have 


^02 ■ 


~ ^i-lj ^ ^i+2 j ^i+lj 
""i ■ ''i-2 ''i+2 ■ 


(D-9) 


(D-10) 


(D-U) 


Then the equation becomes 

(Coi + Co2/2)[5i] - C „3 < >/2 = * C „3 < >/2} (D-12) 

We now express the jump and average values across the shock in terms of the values of (p at the grid points 
for the point for which Uj ^ ^ j > 0 and u^j < 0. We have 




(D-13) 


<‘^lx > 


1 r *^i+l j ~ ‘^ij + ‘ ^i-1 j ~ *^1-2 j 1 

2 \ - X. xj.j - X..2 J 


(D-14) 


WeletCo4 = Co3/4(Xi^.i-Xi)andCo5 = Co3/4(Xi.i - x^.g). Then we get 


(Cqi + Co2/2)(v>ij <pj.j j) Cq4(v>j^j ■ - (pjj) - Cq5(v>j,j j - j) 

= ^l[(Coi - Co 2''2)(‘'« - <.li) j - %) 


(D-15) 


This equation replaces the difference equation in the x sweep. 
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APPENDIX E 

A CONJUGATE GRADIENT ALGORITHM FOR 
ASYMMETRIC INDEFINITE COMPLEX MATRICES 


In this section, we present a conjugate-gradient-type algorithm for unsymmetric complex matrices. The 
algorithm has an asymptotic convergence rate inversely proportional to the square root of the condition number 
of the coefficient matrix and does not assume the coefficient matrix to have any special property other than 
nonsingularity. 

In section E.l, we shall present the principle of the classical conjugate gradient method. In section E.2, we 
shall present our algorithm and related theories. 

El DESCRIPTION OF THE CONJUGATE GRADIENT METHOD 

Suppose we are interested in finding the solution x of the equation 


Ax = b (E-1) 

In the case that A is symmetric (or Hermitian in the case of complex matrices) and positive definite, the k-th step 
of the classical conjugate gradient method finds an approximate solution x^ for equation (E-1), so that the norm of 
the residual is minimum in the space spanned by the vectors { Av^, Av 2 , . . AV|^} , where Vj, V 2 , . . V|^ are 
orthonormal vectors, that is, (v^, Vj) = 0, for i ^ j, and the lengths of the v^’s are equal to one. These orthonormal 
vectors are the backbone of the conjugate gradient method, and is a linear combination of these Vj’s. 

Provided these Vj’s can be computed easily, there is a whole class of algorithms which minimize the norm of 
the residual in the vector space spanned by Av^, Av 2 , . . . Av^ in the k-th step. These algorithms can be considered 
as variants of the classical conjugate gradient method. Among these are the conjugate residual, the modified 
conjugate residual, the variational method, SYMMLQ (Conjugate gradient method with LQ factorization for 
symmetric indefinite system), LSQR (Conjugate gradient method with QR factorization for least squares prob- 
lems, A* A X = A*b, where A is rectangular matrix), etc. Chandra (ref 33) gave a detailed comparison of the ones 
applicable to symmetric matrices. 

For the sake of completeness, we present the classical conjugate gradient algorithm: 

Algorithm 1. Classical Conjugate Gradient 

Step 1: Choose Xq 

Compute ro = b-Axg 
Set po = ro 

ao = (ro,ro)/(Po,Apo) 

Xi = Xq + aoPo 
ri = ro - aoApo 

Step 2: Compute 

ai = (rj,ri)/(pi,Api), 

Xi -I- 1 = Xi + ajPj, 

^i + i = i-i-aiApi 
bi = (ri + i,i-i + i)/(ri,ri) 

Pi + i = Ti + i + biPi 
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Step 3: If ^ 1 is small enough, terminate the iteration; else set i = i + 1 and go to step 2. 

The vectors Pj are orthogonal to each other, although they are not unit vectors. The vectors v^ may be constructed 
by normalizing p^. 

In section E.1.1, we shall discuss orthogonalization. In section E.1.2, we shall discuss solution of equation 
(E-1) based on the orthonormal vectors. 

E.1.1. ORTHOGONALIZATION 

There are two classical orthogonalization processes for solutions of different numerical linear algebra 
problems. 

First, there is the Lanczos’ orthogonalization process for symmetric (or Hermitian) positive definite 
matrices. This process is computationally very efficient and extremely economical because it only requires the 
presence of Vj^ and Vj^.^ to compute v^^ ^ j. If we let Vj^ be the matrix whose columns are v^, V 2 , . . V|^, then the 
Lanczos’ orthogonalization process at the k-th step has performed the following decomposition: 

(Vk*)AVk = Tk (E-2) 


where T|^ is a tridiagonal matrix. 

More general for asymmetric matrices is the Gram-Schmidt orthogonalization process. This procedure 
does not have the computational efficiency of the Lanczos’ process. To find the Vj^ vector, it requires presence of all 
the Vj’s computed before. At the k-th step, the Gram-Schmidt orthogonalization process performs the following 
decomposition: 


(Vk*)AVk = Rk (E-3) 

where is an upper triangular matrix. When the dimension of A is large, we see that the Gram-Schmidt proc- 

ess is computationally infeasible. 

The new algorithm we present here involves finding two sets of orthonormal vectors U and V. Similar to 
the Lanczos’ process, it only requires the presence of v^.i and v^, and u^.i and u^ to compute v^ + 1 and u^ + ^ The 
k-th step of our algorithm performs the following decomposition: 

Uk*AVk = Tk (E-4) 

E J.2 COMPUTATION OF THE APPROXIMATE SOLUTION 

Note that the second equation of Step 2 in Algorithm 1 computes the (i + l)-th approximation of the solution. 
The vector p^ is sometimes called the directional vector. Other approaches proposed by Paige and Saunders (ref 
40), Bunch and Kaufman (ref 41), and Chandra (ref 33) to extend the classical conjugate gradient method to 
symmetric (Hermitian) nondefinite problems have been motivated by the following observation: 
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Let A be a symmetric (Hermitian) positive definite matrix. In the k-th step of Lanczos’ orthogonalization 
process, the following matrix decomposition is performed: 

Vk Wk = T, 

If we choose = ro/||ro || , then x^. in equation (E-5) below is mathematically the same as the approximation 

to X in the k-th iteration of the classical conjugate gradient method: 

Tkyk = Ikollei (E-5.a) 

and 

Xk = xo + V^yk (E-5.b) 

(e^ is the vector with 1 in its first entry and zeroes everywhere else.) 

Of course when A is symmetric (Hermitian) but not positive definite, the above statement is no longer true; 
but one can still solve equation (E-5) to obtain an approximate solution. 

We have found that solving equation (E-5.a) by an LQ factorization (where L is lower triangular and Q is 
orthonormal) is numerically much more stable than generating directional vectors for the approximation of Xj^’s. 

Thus there are two features in our new algorithm that are distinct from the classical conjugate gradient: 

1. A new way of generating orthonormal vectors. 

2. A tridiagonal system by an LQ factorization to obtain an approximation to the solution. 

E1.3. TRIDIAGONALIZATION ALGORITHM FOR UNSYMMETRIC MATRICES (USYMLQ) 

As stated in section E.1.2, the k-th step, k = 1,2, ... of the algorithm to be described in this section involves 
the following decomposition: 


Uk*AVk = Tk 

where Uk and Vk are matrices whose columns are orthonormal and Tk is a tridiagonal matrix. We shall first 
describe our process with which we generate the orthonormal vectors, the Uk’s and the Vk’s: 

Algorithm 2 (Tridiagonal Process) 

(a) Set Uq = 0, 

PiUi = b,YjVi = c 

(b) 

(b.l) p = Auj - YiUj.i 

(b.2) q = A*Uj - PjVj.i 

(b.3) Oj = u*iP i = l, 2, 3, ... 

(b.4) Pi+iUi + i = p-UiUj 

(b.5) Yi + iVj + i = q-a*;Vi 
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where b is the right hand side vector of the matrix equation we want to solve, and Pj and are chosen so that the 
vectors and Vj will have norms equal to 1. We shall discuss the choice of c later in this section. For the moment, 
let us set c = b. The process terminates when P^ or y^ equal zero. 

In the rest of this section, we shall give an overall picture which relates the k-th approximate solution Xj^ to 
the k-th step of Algorithm 2. If we substitute (b.l) and (b.2) into equations (b.4) and (b.5) respectively and let 
and V,^ be nxk matrices whose columns are respectively the u^'s and the v/s obtained from the first k steps of 
algorithm 2, and if we define a tridiagonal matrix as 


Tu = 


«iY2 
P 2 «2 Ya 

Ps 


••••Yk 
3k «k 


where Pj, and y^ are defined for each i in equation (b), then the first k steps of Algorithm 2 can be written in the 

following matrix equations: 


AVk “ Uk^k + Pk+i^k+A^ (E-6.a) 

A*Uk = VkTk* + yk.iVk.iekT (E-6.b) 

where ek is a vector of length of n with 1 at the k-th entry and zero everywhere else. Multiplying (E-6.a) by an 
arbitrary k- vector, yk, whose i-th element is ri^, we obtain 

^"^kyk = UkTkyk + Pk + iUkOk'^yk 

Since b = UkCPiCi) by definition, then if yk and Xk are defined by the following equations 

Tkyk = PiCi (E-7.a) 

Xk = Vkyk (E-7.b) 

then we shall have 


Axk = b -I- ilkPk + iUk+i 

Hence Xk may be taken as an approximation to x, and will solve the original system if rikPk + 1 is negligibly small. 
The above arguments are not complete, but they provide some motivation for defining the sequence of vectors Xk 
according to equations (E-7.a) and (E-7.b). 

If Xk is as defined in equation (E-7), then the residual rk will be 

rk = b - AXk (E-8.a) 


“ Uk + iCPxe^-Skyk) 


(E-8.b) 



where 



If we define = P^e^ - U^ + iS^yk, then = U^ + itk, and Hr^ll = |ltk||, assuming 11^ + 1 has been computed by 

exact arithmetic. Hence it is natural to solve the least square problem 

min IIPiei-TkykIl (E-9) 

yk 

which is the basis for USYMLQ. Computationally it is advantageous to solve this problem by the standard LQ 
factorization (see Paige and Saunders ref. 42). 

In this section, we have given a brief discussion of the motivation behind the algorithm USYMLQ. We shall 
go into greater details concerning the theoretical and computational aspect of the algorithm. 

E.2 THE VECTORS OF ALGORITHM 2 

In section E.2.1, we shall discuss the properties of the vectors u^. and v^ generated by Algorithm 2; in sec- 
tion E.2.2, we shall discuss the possible choices for the initial vectors % and v^; in section E.2.3, we shall discuss 
in detail the LQ factorization. The arguments used for the proofs of the theorems in these sections are based on 
basic linear algebra tools such as linear independence, minimum polynomial, etc. These sections are included 
for the sake of completeness. 

E.2.1 PROPERTIES OF THE VECTORS AND V^ 

This section describes the theory governing the vectors u^ and v^ generated by Algorithm 2. 

We define the symbol < x^, X 2 , . . ., x^^ > as the vector space spanned by the vectors x^, X 2 , . . x^^. The following 
theorem summarizes the properties of and v^ which are generated by Algorithm 2. It is proved by induction. 

Theorem 1. The following four statements are true: 

1) Uk*Uk = I 
V,*Vk = I 

2) U 2 k is a linear combination of vectors in the vector space 
U 2 k = < b, AA*b, . . ., (AA*)k-ib, Ac, AA*c, . . ., (AA*)k-iAc > 

U 2 k-i is a linear combination of vectors in the vector space 
U 2 k + i = <b,AA*b,...,(AA*)kb,Ac,AA*c,...,(AA*)i^-iAc> 

Vgk is a linear combination of vectors in the vector space 

= < c, A*Ac, . . (A*A)k-ic, A*b, A*Ab, . . (A*A)k-iA*b > 

V 2 k + 1 is a linear combination of vectors in the vector space 

^2k + l 

= <c, A* Ac, .... (A*A)kc, A*b, A*Ab, . . (A*A)kiA*b> 
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3) Let A = UDV* be the singular value decomposition of A; that is, tJ and V are orthonormal matrices, and D is 
a positive and real diagonal matrix: D = diagtdp . . ., d„) and let 


b= j^jhiUi, 

kiVi, 


i = l 


where the hj and the are scalars, is the i-th column of the matrix U and Vj is the i-th column of the matrix 
V. Ifb = c, then the following is true: 

C/2k = < i^di^h^U,,..., J^di2(k-l>h^Ui, i^d^k^Ui,..., i^di2<k-l).^k.a.> 

f^2k + l = < ilihiUj, Jjdi2hiUi,..., .|jdi2khjUi, J^dikjUi,..., J^di2<kl) + lkiUi ) 

^2k = < Jjdi2kiVi,..„ i^di2(k-l)k.^j, J^dihjVi i^di2<k-l) + lhiVi> 

^2k + i = < J^di^kiVi.-”. jljdi^'^kiVi, J^djhjVi, ..., J di2<k-i>+ihiVi > 

Also if X|^ is generated by (E-7.a) and (E-7.b), then the k-th residual rj^ = b - Ax^ is in t/k+i- We use the 
symbols and Vj^ (without the hats) to denote the orthonormal vectors generated by Algorithm 2, and Uj^ and 

(with the hats) to denote the singular vectors of the matrix A. 

4) If is computed by equation (E-7), and the residual vector is as defined by equation (E-8), then r,^ 
implies n^ + i^O and Vj^ ^ ^ 0. 


Note that parts 1) and 2) of the Theorem 1 described the “Lanczos-like” characteristics of Algorithm 2. Part 

3) of Theorem 1 implies that the asymptotic convergence rate of the conjugate-gradient-type method based on 
Algorithm 2 is dependent on the square root of the condition number of A rather than the condition number 
itself, as in the case of solving the normal equation A*Ax = A*b by the classical conjugate gradient method. Part 

4) takes care of the case when Algorithm 2 terminates, that is, when Uj^ or v,^ are zero vectors; this part of the 
theorem asserts that in the case when the algorithm terminates, we have arrived at the solution. 

We called the procedure of computing the Uj^’s and the Vj^’s by Algorithm 2 and the Xj^’s by equations (E-7.a) 
and(E-7.b)USYMLQ. 

Proof of Theorem 1 

1) We prove the following statements by induction: 


(Ui,Uj) = 0 fori<j (E-lO.a) 

(vj,Vj) = 0 fori<j (E-lO.b) 

(Avi,Uj) = 0 for i<j -I- 1 (E-lO.c) 

(A*Ui,Vj) = 0 for i<j - 1 (E-lO.d) 


When i = 1, the statement is true by construction. Suppose the statement is true for i = k; we show that it is true 
for i = k H- 1. By Algorithm 2, 


Pk + 1% + 1 = Avk - YkUk-i - «kUk 


59 



Forming inner products with Uj on both side of this equation, we obtain 

Pk + i(Uk + 1. Uj) = (Avk, Uj) - Yk(Uk.i,Uj) - ak(Uk,Uj) 

When j < k-1, by the induction hypothesis regarding (E-lO.a) and (E-lO.c), we see that (u^ + i,Uj) = 0. When j = k, the 
definition of and the induction hypothesis regarding (E-lO.a), (u^ + i,Uj) = 0. As for j = k-1, we note that the last 

line of part (b) of Algorithm 2 yields: 


A*Uk.i = YfcVk + PkVk.2 + c(k*Vk.i 


which again implies: 


(Vk,A*Uk.i) = Yk(Vk,Vk) + Pk(Vk,Vk.2) + ak*(Vk»Vk.i)- 

By the induction h 3 rpothesis regarding (E-lO.b), the equality we just derived yields: 

(Vk,A*Uk.i) = Yk(Vk,Vk) = Yk 

The properties ofthe inner product imply (Vk,A*u,^-i) = (Avk,Uk,i).Thus(Uk + i,Uj) = 0whenj = k-L Therefore we 

have proved (E-lO.a). The proof of (E40.b) is similar and therefore will not be repeated here. By construction 
(uj^,Uj^) = land(Vk,Vi^) = Iforallk. Therefore the proof ofpartl)ofTheoremlwillbecompletedifweprove(E-lO.c) 
and (E-lO.d) above. 

Note that Algorithm 2 implies 


AVk + i = Pk + 2Uk + 2 + Yk + lUk + Otk + iUk + x, 

so 

(Avk + i.Uj) = Pk a(Uk + 2-Uj) + Yk + i(Uk,Uj) -H Uk + i(Uk + i.Uj) 

which equals zero when j<k. The proof of (E-lO.d) is similar and therefore is omitted here. Thus we have com- 
pleted the proof of part 1) of Theorem 1. 

2) The proof of part 2) of Theorem 1 is also by induction. When k = 0, the statement in 2) is true by construction. 
Now we assume the same statement is true for k and prove that it is true for k + 1. Note that for any non-negative 
integer j, 

Uj 

Vj 

A*Vj CV.^, 

AVj 

Thus Pk + lUk + 1 = Avk - OkUk - Yk^k-i s 14 + 1 - Similarly, we can establish the proof regarding Vk + j. 

3) The singular decomposition A = UDV* implies that 

AV = DUandA*U = DV 
which is equivalent to 
Avj = dj Uj and A*Uj = d^Vj for j = 1, 2. . 
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Since 


and 


Then 


and 


Also 


and 


b = £ h:li: = I k:V:, 

j=l J J j = l J J 

Ab = IjkjAvj = J^kjd^Uj, 
A*b = J^hjA*Uj = .i hjdjVj 


A*Ab= J^kjdjA*Uj= J^kjdj2vj 
AA*b = J^hjdjAdjVj = J^hjdj2uj 
(A*A)'^b= L ki(L2kVi 

j=l J J J 


(AA*)kb = L^hjdj2kuj 

Also note that 

(A*A)*'A*b = A*(AA*)><b = E h:i2kA*Ui = E h:cL2k + iv. 

j=i J J J j=i J J j 


and 


(AA*)kAb = A(A*A)kb = J^hjdj2kAvj = ,E hjdj2k + i^i. 

Substituting all these equalities into part 2) of Theorem 1 will yield part 3). 

4) To prove part 4), we are going to show that if either u ^ + 1 or + j, the Gc + D-th orthonormal vector generated by 

Algorithm 2 is zero; then r^ = 0. Let and be nxk matrices whose columns are respectively the Uj’s and the 

Vj’s in the first k steps of Algorithm 2, and Tj^ is the tridiagonal matrix which was defined above equation (E-6.a) 
and (E-6.b). If we multiply both sides of equation (E-6.a) on the right by the matrix U|j, then 

U,*AV, = T, 


Let r^ be defined in equation (E-8); then 


UkX = Uk*(b-Axk). 
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Since x,, is as defined by equation (E-7) and b = PiUi^e^ by Algorithm 2, the above equality can be rewritten as 


UkX = Uk*(UkPiei-AVkyk) 

= Piei-Uk*AVkyk 
= Piei-Tkyk 

which is equal to 0 by construction. Thus (uj,rk) = 0 for j < k. However, equation (E-8) implies that 

Uk + 1% = Uk.,i*(Uk.,i*(3iei-Skyk)) = PA-Skyk = tk- 

If U|,^.i = 0, then tj, = 0. However as we have pointed out in a remark following equation (E-7), 
ll^kll = II tj^ll- Therefore tj^ = 0 implies rj^ = 0. 

We shall now consider the case v^^ ^ ^ = 0. If we define Xj^ and yj^ by 


Tk*yk = Pa (E-n,a) 

Xk = Ukyk (E-E.b) 

then x,^ will be an approximation to the solution of A*x = b. Using an argument similar to the one used for the 
caseU|^ + i = 0,wecanshowthatifV|^^i = 0,thenfi^ = b-A*Xj^ = 0. However, we are not interested inf and we 
really want to show that rj^ = 0. 

Note that f = 0 means the solution x for A*x = b is in and rj^ = 0 means that the solution x for 
Ax = b is in V^. We shall show that x in implies x in Let us first write these solution vectors in terms of 

the singular vectors and singular values of A. Recall that 

b = z hjUi = Uh 

I kw: = Vk 

i = l * ^ 

where h is a vector whose entries are h^, h 2 , . . and k is a vector whose entries are kj, k 2 , . . .. Note that 

X = A-ib = VD itJ*b = VD iU*Uh = VD % = £ d-ihiVi. 

i = l * 

Similarly, 


X = (A*)-ib = J^di-ikiUi 
For simplicity, let us assume k even and replace k by 2k. x £^ 2 ^ means 


X = E dfikiUi 

i=i ^ 

= il/ + *'ionjdi2j+iki)Ui 

Since the Uj’s are linearly independent, this means that for i = 1, 2 

dj ^ki = ( ^E^m^dj^ihi + n^dj^j+ikj) 
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Multiply both sides of the above equation by dj, and move the left hand side to the right hand side; we obtain 


0 = ( L mjdj2j + ihi + E njA^Jk;) 
j=o J ' ‘ j=o-) ' * 

where n„ = -1, and nj = nj + j for j = 1, 2 

Multiply both the above equations on the right by djUj and add all the equations together: 


0 = 



The vector on the right hand side of the above equation is in f/gk + 1- This equality implies that the basis vectors of 
^2k + 1 defined in part 3) of Theorem 1 are linear dependent. Therefore, we can express the first basic vector 

n 

L h:U: as a linear combination of the other basis vectors: 

i=l ‘ ' 


Z hiUi = 

i = l * ‘ 



The linear independence of the vectors implies that for each i = 1, 2, . . 


hi = ( \^m:di2jhi + \^nid:2j + lk:X 
* j=l J * » j=o J ' ' 

Multiply both sides of the equation by and sum each side to obtain 


,E dfiVih; = 



Note that the left hand side is the solution vector for Ax = b, and the right hand side is a vector in Vgk- The 
conjugate gradient algorithm can be considered as choosing Xak in so that the residual norm is the least. 

Therefore, we conclude that r2k = Owhenvgk+i = 0. 


The proof for k odd is similar and therefore is omitted here. 

The proof of part 4) of Theorem 1 yields the following corollary which is useful in defining a stopping crite- 
ria for our algorithm: 

Corollary: ||rkl| = HS^yk - Ple||, where r^ and are as defined in equation (E-8). 


E.2,2 CHOICE OF THE INITIAL VECTORS Uj AND Vj 

Choosing u^ siich that P^Ui = b is reasonable. Such a choice is made in the standard CG-type algorithms 
like SYMMLQ and LSQR. However, we know of no works in the literature to provide us with guidelines or 
motivation regarding the choice of v^. The proof of part 4) of Theorem 1 shows that a solution procedure for the 
equation A*x = c is “implicitly” embedded in Algorithm 2, combined with equations (E-7.a) and (E-7.b), and this 
implicit procedure occurs simultaneously with our explicit procedure for Ax = b. The proof of part 4) of the the- 
orem also implies that we want to choose c such that we do not solve A*x = b before we solve Ax = b. This is why 
we cannot choose c aA*b, because the solution vector, which happens, for A*x = c, to be b, lies in U^, A*x = c is 
solved (implicitly) in 1 step, Algorithm 2 terminates and we are not even close to the solution for Ax = b. In the 
proof of part 4) of Theorem 1, we also showed that if we choose c = b, then A*x = c cannot be solved (implicitly) 
before Ax = b is solved. In other words the choice of c = b guarantees that our algorithm will arrive at the solu- 
tion. However, this is not the optimal choice. The reason for this is given in the following theorem. 
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Theorem 2. If the matrix A has m distinct singular values, then the number of iterations USYMLQ takes to 
converge is bounded by min(2*m,n), while the number of iterations required by the classical CG applied to the 
normal equation is bounded by min(m,n). 

Proof of Theorem 2 


From part 3) of Theorem 1, we see that 


V2m = <c, A*Ac, . . (A*A)^-ic, A*b, A*Ab, . . (A*A>”'^A*b> 

= < J^kiVi,di2kiVi, J^di2(m-i)k.v.^ 

. djhjVi , . . . , d. 2(m-i) + ihj vkj > 

which is obtained after 2m steps, contains the solution vector. As for applying the classical CG to the normal 
equations, we see that the solution vector is in 

= < JiVi. .IjdjZkiVi,..., 

Basis vectors for are generated within k steps of the classical CG applied to the normal equation. Thus the 
theorem is proved. 


We note that the solution vector is contained in the vector space: 


= < I k:W 


i = l 


.E^dikiVj, 




n II ^ 

Note b = L hiUj, for some h/s. If we could choose c such that c = ^he number of 

iterations required by USYMLQ will be bounded by min(m,n). With an argument similar to the proof used for 
part 4) of theorem 4, we can show that we will not solve A*x = c before we solve Ax = b. However, in practice, we 
do not want to compute the singular values and singular vectors of the matrix A. Although we can express c in 
terms of A, A*, and b, namely, c = A KAA *)‘^"b, we know of no computationally efficient way of computing this 
ideal c. 


The above discussion regarding this ideal c proves the following theorem: 


Theorem 3 


If c = A*HAA *)‘^^b, and if Uj^ and Vj^ are as defined in Algorithm 2, and Xi, is as defined in equations (E-7.a) and 
(E-7.b), then r^^O implies u^^O, and also the number of iterations for this algorithm to converge is 

bounded by min(m, n) where m is the number of distinct singular values of A, and n is the order of the matrix A. 
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E.2.3 THE LQ FACTORIZATION 

In this section, we describe the solution of equations (E-7.a) and (E-7.b) in detail. As we remarked before, it 
is computationally advantageous to solve this system by the standard LQ factorization on T^. The LQ factoriza- 
tion on Tj^ takes the form: 


SkQk - 



’ Qk - Lk - 


Pi 

S1P2 
^1^2 P 3 
^283 - 

^3 - 

Pk-1 

8k-lPk 

^kA 


(E-12) 


where Qj^* = Q 2 ,i-Q 3,2 • • • Qk,k-i a product of plane rotation designed to eliminate the superdiagonal elements 
Yi, Y 2 , Therefore, the left hand side of equation (E-7.a) can be written as 

TkYk = Lj^Qk “ Pi^i 

We note that has no elements in common with yj^.^. We are actually not interested in y^., but rather it is as 
defined in equation (E-7.b) that we are after. We want to organize the computation of Xj^ so that only the most 
recent iterates need to be saved. The scheme we discuss below avoids the explicit computation of yj^. We let = 

Qkyk> - VkQk* Thus x,^ = Wj^Zj^. We will show that z^ and Wj^ can be computed by simple recursions. 

If we let = Qk, yk, then the first Gc-2) entries of Zy. are the same as those of z^,^. That is, if we write 
Zk = (zi, Z 2 , . . Zk.i, ZkF, then 


^k-l “ ("^k -3 ^k -3 " 8^.2 2 ^. 2 )/ Pk-I (E-13.a) 

^k “ ("^k-2 ^k-2 ‘ 8k-i Zk-iV Pk (E-13.b) 

Note Zk will be replaced by Zk in the next iteration. 

Now consider the plane rotation Qk,k i- Qk,k-i operates on the (k-l)-th and k-th row of Tk to eliminate Yk- This 
gives the following simple recursion: 



(E-14) 


where and = p^, and the scalar Ck and Sk are the nontrivial elements of Qk k-i* The quantity pk and 6k 

are to be replaced by Pk and 6k in the next iteration. Consider the nxk matrix 

Wk = VkQk* 
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Note that the first (k-2) columns of Wj^ are the same as those of That is, if we write 


Wk = (Wi,W2,...,Wk.i,Wk) 

where the w/s are the columns of Wj^, then 


Wk.i = Ck* Wk.i + SkVk (E-15.a) 

^ = SkWk.i-CkVk (E45.b) 

Note that will be replaced by in the next iteration. Therefore, if we first substitute = Q^y^, then 
into equation (E-7.b), we get 


If we let Xq = 0, 


Xk = VkYk = VkQkX = WkZk 


^k-1 ^k-2 ^k-l^k-l» 


(E-16) 


then 


Xk = Xk-i + ^^kZk 


(E-17) 


Note that w,^,i, and x together with z are computed in the k-th iteration. Also note that is to be 
replaced in the next iteration. We could have computed according to equation (E-17), but it is not necessary to 

do so until the last iteration. From the above discussion, we see that only the last two iterates need to be saved. 

Now we consider the stopping criteria of our algorithm. Two stopping criteria are implied by part 4) of 
Theorem 1 as well as Corollary 1: 

a) If pj^ or Yj^ are small enough, rj^ will be small enough. 

b) However it is possible for r^^ to be small when the norms of either u^ or v^ are not small. Corollary 1 permits us 
to compute the norm of the k-th residual rj^ explicitly: 

llfkll = IISk ■ PaII = ISkZk.i - Ck Zkl*Pk + 1 (E-18) 

In summary, at the k-th iterations, the following quantities are computed: 

i) The vectors u^ and Vj^ according to Algorithm 2 

ii) Ck and s^, which are the nontrivial elements of the plane rotation Qk,k-i defined in equations (E-14) and 
(E-12) 

iii) The scalars p^.i, Pk which are the nonzeroes of the last two columns of defined in 

equation (E-12) 

(iv) Zk_i as defined in equation (E-13) 

(v) Wk,i as defined in equation (E-15) 

(vi) Xk_i as defined in equation (E-16) 

(vii) lirkll according to equation (E-18). 
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Note that if the algorithm terminates, we need to compute according to equation (E-17) before returning to the 

calling program. 

Algorithm 3 USYMLQ 

1. (Initialize.) 

PiUi = b,Vi = Ui.Wj = VpXo = 0,Yi = Pi, Pi = ai,6i = Pi 

2. For i = 1, 2, . . ., repeat steps 3 to 7. 

3. (Continue the tridiagonalization.) 

p = AUi-YiUi.i 

q = A*U; - Pi*Vi.i 
= u*p 

Pi + lUi + 1 = p-«iUi 
Yi + iVj + 1 = q-CiVi 

4. (Construct and apply the plane rotation.) 

diagonalelement ofLj inrow(i-l): pj.i = + y. 2 )/ 

plane rotation: Cj = p^.j/ pj.^ 

Si = Yi/Pi-1 

subdiagonal element ofLj in row i:5i = + YiSi 

diagonal element of in row i: 8^ = S^.^s^ - y^Ci 

5. (Update z, x, and w.) 

^i-l “ Pi-1 

Zi = (-^i-2Zi-2-Si.iZi,i)/Pi 

Xi-l = Xj,2 + Zi.i(Ci*Wi_2 + SiVj) 

Wj.i = Wi.iSj-ViC; 

6. (Test for convergence.) 

If one of the following holds 

(i) Pj is small enough 

(ii) Yi is small enough 

then Wi = Wj.^Si - v^Si 

Xi = Xi + ZiWi 

EXIT 

else go to step 3. 
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E.3 NUMERICAL EXPERIMENTS 


Our new algorithm, USYMLQ (UnSYMmetric LQ factorization) looks very much like Paige and Saunders’ 
LSQR (Least Square QR factorization) (see ref 40). The theoretical results in sections E.l and E.2 show that they 
are indeed different algorithms. In this session, we present numerical experiments which demonstrate that they 
have different convergence properties. 

E.31 GENERATION OF PROBLEMS 

We generated three test problems to compare the performance of USYMLQ with LSQR. The test problems 
have coefficient matrices of the form: 


A = PDQ 

where P and Q are Householder transformations and D is a complex diagonal matrix. Since P and Q are 
orthonormal matrices, the condition number of A, say, K, is defined as follows: 


where d^^^^ and are respectively the largest and smallest absolute values of the diagonal entries of D. More- 

over, the singular values of the matrix A are just the absolute values of the diagonal entries of D. 

We set the i-th entry of the true solution, x, for our test problems to be (n-i) where n is the order of the matrix. 
The right hand sides, b, are obtained by multiplying the true solution with the coefficient matrix: b = Ax. 

E.3.2 CONCLUSION FROM THE NUMERICAL EXPERIMENTS 

The numerical results of our experiment are illustrated by the graphs in figures 58, 59, and 60. The x-axis 
of the figures is the iteration number; the y-axis is the residual norms. The solid lines represent the convergence 
history of USYMLQ, and the dashed line represents the convergence history of LSQR. 
N stands for the order of test matrix, K is the condition number, and S is the number of distinct singular values. 

Test problem 1 illustrates that USYMLQ can handle an ill-conditioned problem with distinct 
singular values quite satisfactorily whilst LSQR seems to have trouble. The result of test problem 
2 is predicted by Theorem 2; when there are multiple singular values, LSQR converges better. Test problem 3 
illustrates the performance of the algorithms for large, well-conditioned matrices. In this case, the convergence 
rates approach the asymptotic convergence rates. 

E.4 THE INCOMPLETE FACTORIZATION 

There are many variants of the incomplete factorization of a matrix. The one that we have been using is a 
product of a unit lower and an upper triangular matrix: M = LU, and M satisfies the following properties: 

1. Ifa^ j ^ 0,thenmi j = a^ j. 
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2. L has the same sparsity structure as the lower triangular matrix of A, and U has the same sparsity structure 
as the upper triangular matrix of A, 

where A is the matrix obtained from the finite difference equation of the small disturbance transonic flow poten- 
tial equation, without the wake terms. In other words, in the two-dimensional case, A is a 5-diagonal matrix. We 
obtain L and U by starting with the regular Gaussian elimination, and ignoring all the “fill-ins.” 

As an example in the two-dimensional case, when IMAX = 5 and JMAX = 4, we may write M = LU as in figure 
61 with the coefficients defined by Algorithm 4 below: 

If we write A = (a^ j), then Algorithm 4 describes L and U for the harmonically oscillating flat plate problems: 
Algorithm 4 (Incomplete LU factorization) 

^1 = ^ 1,1 

Fori = l,2,....(IMAX-2)*(JMAX-2) 

For i = 2,3,..., JMAX-2 
di = aj^i-liXi 

Fori = JMAX-1, JMAX,...aMAX-2)*(JMAX-2) 
li = 

Pi “ ^i,i-JMAX-2/di.i 

= aj j - IjUj.i - Piqi-jMAX-2 

Then R = (r^ j) = A - M has one subdiagonal and one superdiagonal of nonzeroes, that is, R is of the form shown 
in figure 62. 

Thus if A is the matrix of the finite difference equation of the small disturbance transonic flow potential equa- 
tion, then 


A-M = R + W 

where W is a matrix with eight nonzero columns contributed by the wake terms. 
We apply the conjugate gradient method to an equation of the form; 

(I + M-KR + W) )x = M*ib 
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Figure 1.— Comparison of ADI With Kernel Function, Pressure Results, Zero Thickness Airfoil, 
M=0.9,k=0.3 
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Figure 2. —Comparison of ADI With Kernel Function, Pressure Results, Zero Thickness Airfoil, 
M=0.9,k=0.45 




ACp 

real 



- 1.0 - 0.6 - 0.2 0.2 0.6 1.0 


O ADI. At = 0.6 

D Compressible 
kernel 


x-axis 


imaginary 



- 1.0 - 0.6 - 0.2 0.2 0.6 1.0 
x-axis 


O ADI, At =0.6 

□ Compressible 
kernel 


Figure 3. — Comparison of ADI With Kernei Function, Pressure Resuits, Zero Thickness Airfoil, 
M=0.9,k^0.6 
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Figure 4.— Comparison of ADI With Direct Solution of Converged ADI Equation, Pressure Results, 
Zero Thickness Airfoil, M=0.9,k=0.6 
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Figure 5.— Comparison of ADI Using Second-Order 6x Representation With Kernel Function, Pressure 
Results, Zero Thickness Airfoil, M-0.9,k=0.6 
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Figure 6.— Comparison of ADi With Kernei Function, Pressure Results for Two Values of Lt, Zero 
Thickness Airfoil, M =0.9, k -0.6 
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Figure 8.— Comparison of ADI With 0PTRAN2, Pressure Results, NACA 64A010 Airfoil, M=0.85, 
k=0.15 (Concluded) 
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FigureQ.— Comparison of ADI With 0PTRAN2, Pressure Resuits, NACA64A010 Airfoil, M=0.86, 
k= 0.4 (Concluded) 
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Figure 10. — Comparison of ADI With OPTRAN2, Pressure Results for Two 6x Representations, 
NACA 64A010 Airfoil, M = 0.85, k=0.25 
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Figure 1 1 .—Comparlsor) of ADI With OPTRAN2, Pressure Results for Two dx Representations, 
NACA 64A010 Airfoil, M=0.85, k^OA 
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Figure 12. — Comparison of ADI With OPTRAN2, Potential Results for Two 6x Representations, 
NACA 64A010 Airfoil, M=0.85, k=0.25 
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Figure 14.— Comparison of Potentiai Resuits From the ADI Using Shock Boundary Conditions With 
OPTRAN2 Using a Shock Point Operator, NACA 64A010 Airfoil, M=0.85,k- 0.25 
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Figure 15.— Comparison of Pressure Results From ADI Using Shock Boundary Conditions With 
OPTRAN2 Using a Shock Point Operator, NACA 64A010 Airfoil. M=0.85,k=0.25 
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Figure 16.— Comparison of OPTRAD3 and OPTRAN3 WithRHOiV, Pressure Distributions, Aspect 
Ratio 3 Rectanguiar Wing of Zero Thickness, M=0.9, k=0.1. Root Chord 
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Figure 1 7.— Comparison of OPTRAD3 and OPTRAN3 With RHOIV, Pressure Distributions, Aspect 
Ratio 3 Rectangular Wing of Zero Thickness, M=0.9, k=0.1,r\ = 0.46 
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Figure 18.— Comparison of OPTRAD3 and OPTRAN3 With RHOiV, Pressure Distributions, Aspect 
Ratb3RectanguiarWingofZero Thickness, M^O.9, k-0.1,r] = 0.94 
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Figure 1 9. — Comparison ofXTRANSS With Experimental Result, Steady Pressure Distributions, 
Aspect Ratio 3 Rectangular Wing With Circular Arc Airfoil, M=0.9 
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Figure 19.— Comparison ofXTRANSS With Experimental Result, Steady Pressure Distributions, 
Aspect Ratio 3 Rectangular Wing With Circular Arc Airfoil, M=0.9 (Concluded) 
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Figure 20.— The Normalized Bending Mode for the Aspect Ratio 3 Rectangular Wing With Circular Arc Airfoil 
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Figure 21. —Comparison of OPTFiADS With Experimental Result, Pressure Amplitude Distributions, 
Aspect Ratio 3 Rectangular Wing With Circular Arc Airfoil, M=0.9,k=0.13 
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Figure 21. — Comparison of OPTRAD3 With Experimentai Resuit, Pressure Ampiitude Distributions, 
Aspect Ratio 3 Rectanguiar Wing With Circular Arc Airfoil, M= 0.9, k=0.13 
(Concluded) 
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Figure 22.— Comparison of OPTRAD3 With Experimentai Resuit, Pressure Phase-Angie 

Distributions, Aspect Ratio 3 Rectanguiar Wing With Circular Arc Airfoii, M=0.9, k=0. 13 
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Figure 22. — Comparison of OPTRAD3 With Experimentai Resuit, Pressure Phase-Angle 

Distributions, Aspect Ratio 3 Rectangular Wing With Circular Arc Airfoil, M=0.9, 
k=0.13 (Concluded) 
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Figure 23.— Comparison of OPTRAD3 With Two Experimental Results, Pressure Amplitude 

Distributions, Aspect Ratio 3 Rectangular Wing With Circular Arc Airfoil, M=0.9,k-0.13 
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Figure 23.— Comparison of OPTRAD3 With Two Experimentai Resuits, Pressure Ampiitude 
Distributions, Aspect Ratio 3 Rectanguiar Wing With Circular Arc Airfoil, M=0.9, 
k= 0.1 3 (Concluded) 
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Figure 24.— Comparison of OPTRAD3 With Two Experimentai Resuits, Pressure Phase-Angie 
Distributions, Aspect Ratio 3 Rectanguiar Wing With Circular Arc Airfoil, M=0.9, 
If =0.73 
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Figure 24.— Comparison of OPTRAD3 With Two Experimentai Resuits, Pressure Phase-Angie 
Distributions. Aspect Ratio 3 Rectanguiar Wing With Circular Arc Airfoii, M=0.9, 
k=0.13 (Concluded) 
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Figure 25.— Three-Dimensional Representation of the Pressure Distribution for an Aspect Ratio 3 
Rectangular Wing With Zero Thickness Airfoil, M=0.9,k^0.13 
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Figure 26. —Three-Dimensional Representation of the Pressure Distribution for an Aspect Ratio 3 
Rectangular Wing With Circular Arc Airfoil, M =0.9, k= 0.13 
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Figure 27.— Correlation ofXTRANSS With Experimental Result, Steady Pressure Distributions, 
Aspect Ratio 4 Rectangular Wing With Supercritical Airfoil, M = 0.7,ri=0.31 
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Figure 28.— Correlation ofXTRANSS With Experimental Result, Steady Pressure Distributions, 
Aspect Ratio 4 Rectangular Wing With Supercritical Airfoil, M = 0.7, r]- 0.59 
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Figure 29.— Correlation ofXTRANSS With Experimental Result, Steady Pressure Distributions, 
Aspect Ratio 4 Rectangular Wing With Supercritical Airfoil, M = 0.7, n = 0.81 
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Figure 30.— Correlation of XTRAN3S With Experimental Result, Steady Pressure Distributions, 
Aspect Ratio 4 Rectangular Wing With Supercritical Airfoil, M = 0.7, r} = 0.95 
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Figure 31. — Correlation of OPTRAD3 With Experimental Result, Pressure Amplitude Distributions, 
Aspect Ratio 4 Rectangular Wing With Supercritical Airfoil, M-0.7,k=0.178 
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Figure 31 .—Correlation of OPTRAD3 With Experimental Result, Pressure Amplitude Distributions, 
Aspect Ratio 4 Rectangular Wing With Supercritical Airfoil, M = 0.7, k =0.178 
(Concluded) 
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Figure 32.—Comparison of OPTRAD3 With Experimentai Resuit, Pressure Phase-Angie 

Distributions, Aspect Ratio 4 Rectanguiar Wing With Supercritical Airfoii, M=0.7, 
k=0.178 
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Figure 32.— Comparison of 0PTRAD3 With Experimental Result, Pressure Phase-Angle 

Distributions, Aspect Ratio 4 Rectangular Wing With Supercritical Airfoil, M -0.7, 
k= 0.1 78 (Concluded) 
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Figure 33.— Comparison ofOPTRAD3 With XTFtAN3S, Pressure Ampiitude Distributions, Aspect 
Ratio 4 Rectanguiar Wing With Supercriticai Airfoii, M=0.7,k^0.178 
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Figure 33.— Comparison of OPTRAD3 With XTRAN3S, Pressure Ampiitude Distributions, Aspect 
Ratio 4 Rectanguiar Wing With Supercriticai Airfoii, M =0.7, k =0.1 78 (Concluded) 
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Figure 34. — Comparison of OPTRAD3 With XTRAN3S, Pressure Phase-Angie Distributions, Aspect 
Ratio 4 Rectanguiar Wing With Supercriticai Airfoii, M=0.7,k=0.178 
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Figure 34.— Comparison of OPTRAD3 With XTRAN3S, Pressure Phase-Angie Distributions, Aspect 
Ratio 4 Rectanguiar Wing With Supercriticai Airfoii, M=0.7,k=0.178 (Conciuded) 
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Figure 35.— Correlation ofOPTRAD3 With Experimental Result, Pressure Amplitude Distributions, 
Aspect Ratio 4 Rectangular Wing With Superciitlcal Airfoil, M=0.7, k= 0.356 
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Figure 35. — Correlation of OPTRAD3 With Experimental Result, Pressure Amplitude Distributions, 
Aspect Ratio 4 Rectangular Wing With Supercritical Airfoil, M=0.7,k= 0.356 
(Concluded) 
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Figure 36.— Correlation of OPTRAD3 With Experimental Result, Pressure Phase-Angle 

Distributions, Aspect Ratio 4 Rectangular Wing With Supercritical Airfoil, M-0.7, 
k =0.356 
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Figure 36.— Correlation of 0PTRAD3 With Experimental Result, Pressure Phase-Angle 

Distributions, Aspect Ratio 4 Rectangular Wing With Supercritical Airfoil, M =0.7, 
k= 0.356 (Concluded) 
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Figure 37. — Comparison of OPTRAD3 With XTRAN3S, Pressure AmplihJde Distributions, Aspect 
Ratio 4 Rectangular Wing With Supercritical Airfoil, M=0.7,k=:0. 356 


121 


I^CpI 



x-axis 


O OPTRAD3 
— XTRAN3S 


(c) 0.81 semispan 




- 1.0 - 0.6 - 0.2 0.2 0.6 1.0 
x-axis 


o OPTRAD3 
— XTRAN3S 


(d) 0.95 semispan 


Figure 37.— Comparison of OPTRAD3 With XTRAN3S, Pressure Amplitude Distributions, Aspect 
Ratio 4 Rectanguiar Wing With Supercriticai Airfoii, M = 0.7,k^ 0.356 (Concluded) 
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Figure 38.— Comparison of OPTRAD3 With XTRAN3S, Pressure Phase-Angle Distributions, Aspect 
Ratb 4 Rectangular Wing With Superchtical Airfoil, M = 0.7, k-0. 356 
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Figure 38.— Comparison ofOPTRAD3 With XTRAN3S, Pressure Phase-Angle Distributions. Aspect 
Ratto 4 Rectangular Wing With Supercritical Airfoil, M=0.7, k=0.356 (Concluded) 
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Figure 39. — Example of Swept Wing Coordinate System in the Physical Plane for a Clipped Delta 
Planform 
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Figure 40.— Example of Swept Wing Coordinate System in the Physical Plane for an Untapered 
Swept Wing With Modified Root Geometry 
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Figure 41.—OPTRAD3 Pressure Distributions for an Untapered Wing With Three Angles of Sweep, 
M = 0.9, k=0. 13, Root Chord 
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Figure 42.—OPTRAD3 Pressure Distributions for an Untapered Wing With Three Angles of Sweep, 
M = 0.9,k^0.13,^ = 0.51 
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Figure 43. — OPTRAD3 Pressure Distributions for an Untapered Wing With Three Angies of Sweep, 
M=0.9,k=0.13,rj = 0.93 
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Figure 44.— Comparison of OPTRAD3 With RHOIV, Pressure Distributions for a 30-deg Swept, 
Untapered Wing, M = 0.9, k=0.13. Root Chord 
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Figure 45.— Comparison of OPTRAD3 With flHO/V( Pressure Distributions for a 30-deg Swept, 
Untapered Wing, M = 0.9, k-0.13,ri=0.51 
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Figure 46. — Comparison of OPTRAD3 With RHOiV, Pressure Distributions for a 30-deg Swept, 
Untapered Wing, M = 0.9, k= 0.13, ri= 0.93. 
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Figure 4Z— Comparison of OPTRAD3 With RHOIV, Pressure Distributions for a 45-deg Swept, 
Untapered Wing, M=0.9, k=0.13. Root Chord 
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Figure 48. — Comparison of OPTRAD3 With RHOiV, Pressure Distributions for a 45-deg Swept 
Untapered Wing, M = 0.9, k= 0.1 3, 0.51 
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Figure 49. — Comparison of OPTRAD3 With RHOiV, Pressure Distributions for a 45-deg Swept, 
Untapered Wing, M=0.9, k=0.13,n = 0.93 
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Figure 50.— Comparison ofOPTRADS With RHOiV, Pressure Distributions fora45-deg Swept, 
Untapered Wing, M=0.9,k=0.5, Root Chord 
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Figure 51.— Comparison of OPTRAD3 With RHOiV, Pressure Distributions for a 45<leg Swept, 
Untapered Wing, M-0.9, k-0.5, ^ = 0.57 
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Figure 52.— Comparison of OPTRAD3 With RHOIV, Pressure Distributions of a 45-deg Swept, 
Untapered Wing, M = 0.9, k=0.5, ri==0.93 








Figure 58.— Comparisons of Convergence Performance ofUSYMLQ andLSQR, N= 10. K= 1(fi. 
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Figure 61. — Incomplete LU Factorization 
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Figure 62.— Remainder of Incomplete LU Factorization 
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